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INTERPOLATION-BASED  CONDENSATION  MODEL  REDUCTION 
PART  1:  FREQUENCY  WINDOW  REDUCTION  METHOD  APPLICATION 

TO  STRUCTURAL  ACOUSTICS 


1  Introduction 

The  efficient  and  accurate  solution  of  large  degree  of  freedom  (DOF)  problems  remains 
difficult  for  linear,  dynamic  and  steady  state  analyses.  Many  problems  have  become  too  large 
and  complex  to  solve  in  an  economic  fashion,  especially  if,  for  example,  a  large  number  of  fre¬ 
quency  responses  are  required  for  a  structural  acoustic  analysis.  In  addition,  the  problem  of  the 
interaction  of  scales,  which  can  extend  from  the  micro  to  the  macro-scale,  has  not  been  ade¬ 
quately  addressed.  This  is  especially  true  for  heterogeneous  media  such  as  composite  materials 
and  structures.  For  linear  problems,  the  use  of  sub-structuring  methods  [1,2]  can  be  very  effec¬ 
tive.  However,  when  applied  to  acoustic  steady-state  analyses,  sub-structuring  methods  introduce 
further  approximations  and  still  require  the  solution  of  the  reduced  problem  at  each  desired  fre¬ 
quency. 

Overall,  finite  element  method  (FEM)  programs  have  increased  in  capabilities  and  in  us¬ 
age  with  improvements  in  computational  platforms.  The  application  of  FEM  analyses  to  large 
complex  problems  has  significantly  risen  as  the  need  for  more  detailed  and  sophisticated  analy¬ 
ses  and  designs  has  increased.  Computationally  based  analysis  is  utilized  in  many  fields  where 
full  experimental  analysis  is  difficult  to  obtain.  Thus  computer  modeling  is  relied  on  as  there  is 
no  other  method  of  obtaining  and  testing  designs.  The  costs  of  both  experimental  and  computa¬ 
tional  analyses  are  a  major  driving  factor  for  improved  techniques  for  large  complex  analyses. 
FEM  computational  tools  which  could  reduce  the  computational  costs  of  very  large  complex 
models  and  which  provide  enhanced  analysis  capabilities  are  desperately  needed. 

Flippen  [3,4]  has  developed  a  general-purpose  dynamic  condensation  model  reduction 
(CMR)  theory  which  utilizes  projection  operators  on  the  governing  differential  equations.  This 
approach  allows  general  DOF  reduction  in  spatially  discretized  models  of  heterogeneous  mate¬ 
rial  and  is  not  restricted  to  periodic  media  or  other  homogenization  assumption.  The  CMR 
method  is  applicable  to  general  second  order  differential  equations  in  time.  This  methodology  [5] 
was  implemented  into  a  linear  finite  element  environment  and  applied  to  dynamic  analysis  (time 
domain)  of  transient  elasticity  problems. 
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In  reference  [6],  Flippen  extended  this  condensation  model  reduction  (CMR)  theory  to 
include  frequency  domain  analysis.  The  use  of  projection  operators  on  the  governing  differential 
equations  allowed  interpolation  across  selected  frequency  windows  to  be  accomplished.  Her- 
mitian  interpolation  functions  [1],  which  required  derivative  information  on  the  boundaries  of  the 
frequency  window,  were  used.  The  CMR  method  was  applied  to  structural  acoustics  analysis  and 
was  implemented  into  a  one  dimensional  (ID)  finite  volume  environment.  The  method  was  ap¬ 
plied  to  a  complex  ID  composite  structure  and  good  results  were  obtained  across  wide  frequency 
windows  for  the  reduced  DOF  solutions  as  compared  to  the  full  DOF  solution.  The  eigenvalues 
of  the  reduced  system  were  found  to  match  quite  well  with  those  of  the  full  system,  even  for  sig¬ 
nificant  DOF  reduction.  In  this  particular  paper,  a  method  was  developed  (though  not  imple¬ 
mented)  to  account  for  singularities  caused  by  resonances  using  singularity  decomposition  tech¬ 
niques.  In  a  recent  effort,  Flippen[7]  has  extended  this  work  into  a  general  model  reduction  the¬ 
ory,  which  unites  and  generalizes  the  reduced  basis  and  frequency  window  class  of  methods  as 
submethods. 

In  this  present  work,  the  CMR  method  has  been  extended  to  frequency  windowing  di¬ 
rectly  into  a  2D  and  3D  finite  element  method  (FEM)  environment  and  applied  to  structural 
acoustics.  The  approach  is  referred  to  as  the  Frequency  Window  Reduction  (FWR)  method  for 
degree  of  freedom  reduction  of  spatially  discrete  time-transformed  linear  models.  In  this  case 
frequency  windows  are  chosen  in  which  a  discretized  model  can  be  reduced  in  its  degrees  of 
freedom.  The  method  consists  of  discretizing  a  FEM  model  and  selecting  the  “master”  node 
points  (DOFs)  which  are  to  be  retained  in  the  analysis.  The  contribution  of  the  “slave”  node 
points,  which  are  to  be  condensed  out,  is  reformulated  and  is  used  to  generate  an  interpolation  of 
the  functional  dependence  of  the  “slave”  response  across  a  selected  frequency  window.  A  “re¬ 
duced”  system  operator  now  may  be  constructed,  which  consists  of  the  “master”  DOFs,  and  the 
influence  of  the  interpolated  “slave”  response.  The  reduced  system  response  can  then  solved 
across  a  selected  frequency  interpolation  window. 

In  addition  to  the  Hermitian  interpolation  functions  (two-point  interpolation)  of  the 
“slave”  response  first  employed  with  this  method,  standard  Lagrangian  interpolation  functions 
(linear,  quadratic  and  cubic  not  requiring  derivative  calculations)  have  been  included.  These  lat- 
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ter  functions,  which  do  not  require  derivative  information  at  the  evaluation  points,  are  an  impor¬ 
tant  addition  that  are  far  simpler  to  implement,  very  efficient  and  provide  function  evaluation 
within  the  window  of  interest  that  can  increase  accuracy.  Also  for  the  standard  quadratic  and 
cubic  interpolation  schemes,  the  location  of  the  evaluation  points  that  are  within  the  interior  of 
the  frequency  window  can  be  varied  for  a  one  dimensional  (ID)  frequency  window  (a  line  in  fre¬ 
quency  space).  This  is  an  important  capability  that  can  substantially  increase  the  accuracy  within 
selected  regions  within  the  interpolation  window.  The  use  of  evaluation  points  within  the  win¬ 
dow  was  extended  to  include  higher  order  Hermitian  functions  with  one  and  two  interior  inter¬ 
polation  evaluation  points  for  ID  windows.  This  has  led  to  a  further  improvement  in  accuracy. 

Computational  analyses  were  performed  using  a  program  developed  in  Fortran  90  which 
can  read  in  and  solve  third  party  generated  finite  element  models.  One-  and  two-dimensional  re¬ 
sults  from  the  FWR  program  eu:e  presented  and  show  excellent  results  for  the  “reduced”  system 
solutions  as  compare  to  the  full  DOF  solutions.  Accurate  results  are  obtained  for  large  DOF  re¬ 
duction  and  large  frequency  windows  even  with  the  presence  of  numerous  resonances  (eigenval¬ 
ues)  within  the  window,  which  have  not  been  explicitly  accounted  for  in  the  analysis  as  yet. 


2  Frequency  Window  Reduction  Method  (FWR) 

The  generalize  equation  of  ODEs  can  be  written  in  compact  form  as 

Lu  =  f  (2.1) 

where  L  is  an  n  by  n  matrix  of  operators,  u  is  an  n  by  1  column  matrix  (vector)  of  the  system  re¬ 
sponse  and  f  is  an  n  by  1  column  matrix  (vector)  of  the  system  stimuli.  The  value  of  n  represents 
the  total  number  of  degrees  of  freedom  (DOF)  for  the  system.  In  order  to  examine  the  frequency 
domain  the  exponential  form  of  the  Fourier  transform,  F  =  e'®',  is  applied  to  equation  (2.1). 
The  operator  L  is  now  a  function  of  the  continuous  independent  variable  ©,  which  generally  is 

complex.  Note  that  i  -  4^ .  Equation  (2.1)  can  then  be  rewritten  as 
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and  the  response  and  stimuli  terms  are 

(2.3a) 
(2.3b) 


u  =  F(u) 
f  =  F(  f ) 


Applying  P  ,  a  permutation  matrix  operator,  on  the  operator  matrix  L  in  equation  (2.1)  or 
L  in  equation  (2.2)  as  L  =  such  that  L  (or  L)  can  now  be  partitioned  as 


L  = 


f, 

^12 

^^21 


(2.4) 


where  the  index  “1”  represents  the  “master”  or  independent  degrees  of  freedom  and  the  index 
“2”  represents  the  “slave”  or  dependent  degrees  of  freedom. 


The  term  “master”  degree  of  freedom  refers  to  the  independent  DOFs  of  the  reduced  sys¬ 
tem  in  the  FEM  model.  The  “masters”  DOFs  are  typically  chosen  as  DOFs,  which  have  applied 
loads  or  constraints  such  as  fixed  boundary  conditions  in  an  FEM  model.  In  addition  to  loaded 
and  constrained  DOFs,  the  master  DOFs  may  also  be  chosen  as  important  nodes  represented  in 
the  models;  for  example,  the  interfaces  or  boundaries  between  element  groups  of  differing  mate¬ 
rial  properties  or  element  grouping  midpoints.  They  may  also  be  selected  as  DOFs  where  the  re¬ 
sultant  system  response  is  of  interest  to  examine. 

The  term  “slave”  degree  of  freedom  refers  to  the  dependent  DOFs  or  those  DOFs,  which 
are  to  be  condensed  or  reduced  out  of  the  full  FEM  model.  These  are  DOFs,  which  typically  do 
not  have  applied  loads  or  constraints.  They  are  also  DOFs  which  do  not  contribute  significantly 
to  the  overall  system  response  or  are  not  of  major  interest. 


Further  applying  the  permutation  matrix  to  the  system  stimuli  according  to 
(\  0^ 


P  =  P 


,-i 


a  0 


results  in  the  system  stimuli  being  partitioned  as 


(2.5) 
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(2.6) 


f  =  Pf  = 

Here  the  subscripts  “m”  and  “s”  refer  to  the  “master”  and  “slave”  DOFs  respectively.  The  ma¬ 
trix  a  is  defined  such  that  =  a  and  I  is  the  Identity  Matrix. 


The  a  matrix,  the  “slave-master”  response  matrix,  is  an  operator  matrix,  which  is  used  to 

relate  “slave”  DOF  properties  such  as  constraints,  forces  and  applied  loads  to  the  “master” 
DOFs.  Usually  “slaves”  DOFs  are  chosen  such  that  they  do  not  have  constraints  or  loads.  How¬ 
ever,  for  very  large  complex  FEM  models  in  order  to  reduce  the  total  number  of  DOFs  ,  it  may 
desirable  to  also  condense  out  DOFs  which  have  constraints  or  applied  loads.  Hence  the  a  op¬ 
erator  matrix  can  be  used  to  interpolate  from  the  “master”  DOFs  for  example  an  applied  load  or 
boundary  condition  to  the  “slave”  DOFs.  In  this  present  work  it  has  been  assumed  that  the  a 

operator  matrix  is  zero  based  on  the  choice  of  the  “master”  and  “slave”  DOFs,  However  a  non¬ 
zero  a  operator  matrix  as  well  as  completely  arbitrary  constraint  and  loading  conditions  of  the 

DOFs  will  be  the  subject  of  future  work. 


The  condensation  model  reduction  (CMR)  theory  has  shown  that  the  general  DOF  reduc¬ 
tion  of  an  operator  matrix  is 


Ljj  -I-  L22  P  —  ^(Lii  +  Lj2 


(2.7) 


where  L  = 


L  L  I 

_  "  .  is  the  reordered  /partitioned  operator  matrix  and  the  a  matrix,  the  “slave- 

V^21  ^22  V 


master”  response  matrix  is  defined.  The  goal  is  to  solve  for  the  j8  matrix  in  order  to  construct  the 
reduced  operator. 


Rearranging  equation  (2.7) 

{t^-at,2)p  =  aL„-L2,  (2.8) 
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The  P  matrix  can  be  found  directly  as 


p  =  G'H  (2-9) 

where  the  G  matrix  is  the  matrix  inverse 

G  =  L22-aL,2  (2.10) 

and  the  H  matrix  is 

H  =  aL„-L2,  (2.11) 

Having  determined  the  form  of  P,  the  final  “reduced"  operator,  Lq,  is 

Lo  =  Li,  +  Li2^  (2.12) 

and  hence  the  reduced  system  equation  to  be  solved  is 

L„  =  n'f 

with  JJf  =  f„  where  =  (I  O)  P. 


The  column  matrix  (vector)  u„  represents  the  “master”  or  independent  DOF  system  response. 
Finally  the  “slave”  response  can  be  obtained  from  the  “master”  response  by  reconstruction  via 


The  goal  in  the  general  CMR  method  is  to  reduce  the  number  of  degrees  of  freedom  in  order 
to  allow  efficient  computational  solutions  to  a  problem  but  which  maintains  the  accuracy  of  the 
full  system.  The  primary  consideration  is  the  accuracy  of  the  reduced  solution  as  compared  to  the 
full  system  solution.  Hence  there  will  be  trade-offs  among  the  degree  of  reduction,  i.e.  the 
number  of  DOFs  condensed  out  from  the  full  system,  the  computational  efficiency  and  the  accu¬ 
racy  of  the  reduced  solution. 
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3  FWR  Formulation  for  Mechanical  Systems 


For  the  case  of  a  mechanical  system  equation  (2.1)  is  written  as 


M— ^  +  D— -  +  Su  =  f 
dt^  dt 


(3.1) 


where  M,  D  and  S  are  the  mass,  damping  and  stiffness  matrices  respectively,  u  is  the  system  dis¬ 
placements  and  f  is  the  system  loads.  So  the  matrix  of  operators  L  is  defined  as 


L  =  M-^  +  D-^  +  S  (3.2) 

dt^  dt 

The  frequency  domain  in  obtained  from  the  exponential  form  of  the  Fourier  transform, 
F  =  ,  applied  to  equation  (3.2),  so  that  the  operator  L  becomes 


L  =  +  imD  +  S  (3-3) 

The  mechanical  system  operator  L  is  a  function  of  the  continuous  independent  variable  ®,  and 
equation  (3.1)  can  then  be  explicitly  rewritten  as 

L(m)  u(6))  =  f(n))  (3-4) 

A 

and  the  displacement  and  load  terms  are  u  =  F(u)and  f  =  F(  f ). 

The  application  of  the  permutation  matrix  operator  to  L,  equation  2.3  as  L  =  PLP"‘ 
such  that  it  can  be  partitioned  directly,  results  in  the  reordering  and  partitioning  of  the  M,  D,  and 
S  matrices  as 


M,,' 

"Du 

D12'' 

and 

'Sn 

0  ^ 
^12 

,  D  = 

S  = 

^22  > 

<®21 

®22> 

^^21 

^22  y 

There  are  two  approaches  to  obtain  the  “reordered”  global  matrices  and  partitioned  matri¬ 
ces  based  on  the  number  of  “master”  and  “slave”  DOFs  used  in  the  CMR  method.  The  first  is  to 
directly  assemble  the  global  matrices  from  the  local  matrices  in  the  reordered  nodal/DOF  num¬ 
bering.  This  approach  requires  an  FEM  code,  which  computes  and  assembles  the  global  mass. 
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damping  and  stiffness  matrices  directly  from  the  element  mass,  damping  and  stiffness  local  ma¬ 
trices.  For  this,  a  look-up  table  of  the  relationship  between  the  original  node  numbering  and  the 
“master”  or  reordered  node  numbering  is  necessary.  Then  for  a  given  element  and  element  con¬ 
nectivity  in  the  original  numbering,  a  transformation  is  used  simply  to  put  that  local  matrix 
(original  numbering)  into  the  global  matrix  (  reordered  numbering).  The  second  approach  is  to 
reorder  according  to  the  “master”  node  number  previously  assembled  global  matrices  in  the 
original  node  numbering  via  the  look-up  table.  This  second  approach  has  the  advantage  that  ma¬ 
trices  may  be  imported  from  other  codes  in  order  to  perform  the  CMR  method.  (See  example 
Appendix  A) 

Finally  the  system  operator  L  matrix  equation  is  obtained  from  the  operator  matrix  in 
equation  (3.5)  in  partitioned  form  as 


9 

+  ico 

fl>n 

(  s 

^11 

S  ^ 
^12 

=  -ccr 

+ 

^22^ 

ll>2. 

*^^21 

^22  > 

(3.6) 


The  general  CMR  method  can  now  be  applied  to  the  above  equation  in  order  to  reduce 
the  number  of  degrees  of  freedom. 


4  (Generalize  Formulation  of  the  System  Operator,  L(<7) 

Beginning  with  equation  (3.3)  the  operator  matrix,  L, 

L(fi!))  s  +  ifi)D  +  S  (4.1) 

can  be  transformed  to  a  complex  variable,  ^  (Appendix  B) 

so  equation  (4.1)  becomes 

L  s  q^M  +  qD  +  S  (4.2) 

It  is  assumed  that  the  matrix  a  =  f{q)  and  that  the  matrices  M,  D,  and  S  are  not  functions  of 
the  generalize  variable^ (or  frequency). 
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Equation  (4.2)  may  then  be  rewritten  as  explicitly  as  a  function  of  the  generalize  complex 
variable  q  in  the  CMR  general  DOF  reduction  form  as 

[L„(«)-a(«)L,j(,)] /)(«)=  C<(?)L„(,)-4(,)  (4.3) 

Expanding  in  terms  of  the  partitioned  mass,  damping  and  stiffness  matrices  and  keeping  a  (q)  as 
a  function  of  q  equation  (4.3)  becomes 

[(9^M22  +  9D22  +  S22)-o^(^)(9^M,2  +  9Dj2  +  8,2)]  ^(^)= 

+  ^D„  +  S„)  -  (9^M2,  +  9D2,  +  S21)] 

The  complex  matrix  ^q)  is  then  found  from  equation  (4.3)  as 

fi(q)  =  G-'(9)  H(q)  <4.5) 

with  the  matrices  defined  as 

a(q)  =  t^(q)  -  a(q)  t„(q)  (4.6a) 

and 

H(q)  -  a(q)  LM  -  t„(9)  (4.6b) 

Alternatively  equations  (4.6)  can  be  written  explicitly  in  terms  of  the  partitioned  M,  D, 
and  S  matrices  as 

Giq)  =  (?"M22  +  qD,,  +  S22)  -  a{q)  +  qD,^  +  S,,)  (4.7a) 

H(q)  =  a{q)  (?"M„  +  ^D„  +  S„)  -  +  qD^i  +  S2,)  (4.7b) 

Finally  reduced  system  operator,  equation  (2.12)  Lo(^)  is 

Lo(^)  ~  L,,(^)+ Lj2(^)  P{q)  (4  8) 

=  +  qD,,  +  S„)  +  (^^2  +  «D,2  +  8,2)  P{q) 
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from  which  the  reduced  system  equation  can  be  solved  for  the  response,  u^,  of  the  “master” 
DOFs 


L.  Un,  =  n?  = 

The  “Slave”  response  u^  can  then  be  reconstructed  from  u^  and  p(q)  according  to  the  equation 

u.  =  P{g)  u„  (4.10) 

In  order  to  obtain  the  reduced  system  operator,  equation  (4.8),  Lo(^),  it  is  necessary  to 
determine  P(q)  from  equations  (4.5, 4.6a  and  4.6b).  From  equation  (4.5),  it  is  required  to  obtain 

G~^(q)  at  every  solution  point,  q  ,  in  order  to  determine  jXq)  and  herein  lies  the  work.  The  de¬ 
sire  is  to  minimize  the  work  necessary  in  obtaining  G~\q)  at  every  solution  point  q,  since  the 
inversion  of  G{q)  would  involve  “slave”  DOF  by  “slave”  DOF  sized  matrix  to  invert  at  each 
point.  So  rather  than  determine  G~^(q)  at  evgrv  solution  point,  a  method  has  been  devised  in 
order  to  obtain  an  interpolated  G~\q)  to  be  valid  within  a  specified  frequency  “window”. 

p{q)  =  G-'{q)  H{q)  (4.11) 

Two  methods,  Hermitian  interpolation  (Section  5)  and  Lagrangian  shape  function  inter¬ 
polation  (Section6),  are  proposed  in  order  to  interpolate  G"‘(^)  at  points  within  a  window  based 
on  the  assumption  that  G~^(q)  is  a  relatively  smoothly  varying  function  of  q  within  the  bounds 
of  the  defining  window.  This  assumes  that  any  eigenvalues  are  within  the  defining  window  not 
significantly  effect  the  smooth  interpolation  of  G~^{q) . 

The  interpolation  method  of  G  ‘(q')  allows  the  formulation  and  calculation  of  the  inter¬ 
polation  or  constructor  matrices  prior  to  any  solution.  Hence  the  main  work  involved  in  the  ma¬ 
trix  inversion  and  matrix  interpolation  calculations  can  be  performed  as  pre-processing  steps. 
This  approach  can  significantly  reduce  the  amount  of  computation  required  to  generate  the  re¬ 
duced  system  operator,  Lq,  and  solve  the  reduced  system  equation 
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(4.12) 


since  the  number  of  degrees  of  freedom  have  been  reduced.  The  reduced  computation  required 
for  solution  of  equation  (4.12)  allows  more  frequencies  to  be  examined  within  the  interpolation 
window.  Also  rather  than  L^,  it  is  possible  to  obtain  directly  the  reduced  mass,  damping  and 
stiffness  matrices,  such  that  other  or  different  types  of  solvers  may  be  use  (Appendix  C).  This 
would  allow  the  use  of  third  party  software  to  generate  the  global  matrices  and  to  solve  the 
reduced  system  problem. 


5  Hermitian  Interpolation  Method 

5.1  Bi-Cubic  Hermitian  Interpolation  of  2D  Rectangular  “Window” 

Hermitian  Interpolations  of  G~\q)  require  the  determination  of  G{q)  as  defined  by 
equation  (4.6a)  or  (4.7a),  the  inverse,  G~\q),  and  the  first  and  second  derivatives  of  G~\q)  with 
respect  to  q  at  the  four  points  defining  the  interpolation  window.  A  “rectangular”  window  is 
specified  by  the  four  comers,  qj[qi  where  Wg  and  hg  are  the  window  width  and  height 

respectively  as  defined  in  complex  q-space. 


qi,q2  +  hg  9i  +  Wg,  qi  +  hg 


w 


The  four  vertices  of  the  interpolation  window  are  then 

=  ^1  +  i  <l2 

^2  =  (^1  +  %)  +  i  q^ 

^3  =  (^1  +  %)  +  *  (^2  +  hg) 
^4  =  i  (^2  +  hg) 
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The  Hermitian  interpolations  of  the  first  and  second  derivatives  of  G  \q)  with  respect  to  q  can 
be  found  by  the  following  approximations 


4c~'W] 

A,, 


-G-'M  ^  G-(?) 

dqj 


(5.1.1) 


'(g)] 

dqjdq^ 


G-.(,)  im  o-'(,)  ^  G-'(,)  + 
^  ^ 


(5.1.2) 


where  G  \q)  is  the  matrix  inverse  of  G{q)  as  defined  in  equation  (4.6a). 


From  the  form  of  G{q), 


G(q)  =  L„(«)-a(«)  L„(«)  (5.1.3) 

the  first  and  second  derivatives  of  G(q)  with  respect  to  qj ,  i.e.  q,  and  ^2  oan  be  found  explicitly 

as 


^G(q)  _  f ^Ljq)  _  ,  .  dUM)  _  Mil 


dq 


a{q) 


dq 


dq 


Li2(^) 


(5.1.4) 


d^Gjq) 

dqjdq, 


f 


dq^ 


-  «{q) 


d%2{q)  _  da{q)  dt,^{q)  _ 

dq^  dq  dq 

dccjq)  dt^ijq)  _  d'^(^{q)  f 

dq  dq  dq^ 


dq  dq 
dqj  dq^ 


(5.1.5) 


_  „(  _  2  aL  fa)  _ 


,  2  a(q)  "'y.^  -  2  , 

dq^  dq^  dq  dq 


dq^ 


dq  d^ 
dqj  dq^ 
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Also  from  the  definition  of  ^  =  q^+  iq^  the  derivatives  — —  can  easily  be  found  as 

dqj 


dq 

dq^ 


=  1 


and 


Having  obtained  the  necessary  matrix  terms,  it  is  now  possible  to  formulate  the  matrices  needed 
to  construct  the  interpolated  G  {q) . 


The  real  interpolation  coefficient  matrix  )  is  calculated  from  the  matrices  )  at 
the  comers  of  the  rectangular  window,  qj  as 


B[q)  =  B[q„q^,w^,\)= 


(5.1.6) 


from  which  the  matrix  inverse  5  is  found.  The  Cj  interpolation  coefficient  matrices  are  ob¬ 
tained  at  each  of  the  points  q^  from 


^1 

^1' 

^2  ^1^2  ^2 

0 

1 

2^1 

3qf 

0 

^2  2^1^2 

^qhi 

0 

0 

0 

0 

1 

qx 

ql 

.0 

0 

0 

0 

0 

1  2q, 

ql 

^1^2 

^2  2 
qyqi 

qUl  ql 

^1^2 

^1^2 

0 

3qUl  0 

^2 

2<?,^2 

3^f^2 

2^2 

'^qUi 

2^?!  ^2  2^2 

39, ?2 

3^f^2 

3^f^2 

0 

2^2 

Hq2 

^qUi  0 

3^2' 

6^,^2 

(5.1.7) 


Defining  the  matrix  Uj[q^  at  a  specific  point  q^  as 
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the  full 


«;(^P^2)= 


G  \qvq2) 

dG~\q^,q2) 

dqi 

^qi 

dq^dq^  ) 


(5.1.8) 


)  matrix  at  all  the  interpolation  window  points  q^j  =  1,2,3 ,4  is  defined  as 


m)= 


g;'(«„32)  ' 

dG~\q^,q^ 

dq^ 

d^G-\q„q,) 

dq^dq^ 

G~'{qi  +M>,q^) 
dG~\q^+w,q^) 
dqi 

dG~\q,  +w,q^) 
dq2 

d‘^G~\q^+w,q^) 

dq^dq^ 

G~\q^-\-w,q^+h) 

+w,^2  +  ^) 

dq^ 

^G"‘(^,  +w,^2  +  ^) 

dq^ 

d^G  +  w,^2 

dqidq^ 

G-\q,+w,q^+h) 
d  G“’(^i  +  w,^2  + 
dq, 

^G“‘(^,  +  W,^2  +^) 
^^G  +  w,^2 


(5.1.9) 
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Finally  ,  the  interpolation  constructor  matrices,  can  be  found  as 

W{q,)  =  "  £/(«,)  (5.1.10) 

NOTE:  Equation  5.1.10  may  be  written  explicitly  as 

;=i 

where  i  =  1,16  and  k  =  1,2, 3, 4,  i,e.  the  interpolation  points.  So,  for  example,  the  first  term  would 
be 


Wife)  -  B,,,  a  («,)  +  fi,.,  +  B,.,  +  ,.4 

A. 5  G  (9=)  +  Si.^— 

n-  C-'(J)  ,  B-  1  B-'  +  B-'  + 

B-'  C-'fe)  I  fl-'  I  JT'  I  B-' 

B, .„  G  (q,)  +  +  ^..is  +  A.,. 

where  fij’  are  the  individual  B  matrix  inverse  terms. 


'  It  is  important  to  note  that  M^(^y  )  and  t^(^y  )  are  16  COMPLEX  “Slave”  DOF  by  “Slave” 
DOF  matrices  and  is  a  16  by  16  REAL  matrix.  The  interpolation  constructor  ma¬ 

trices  MUST  be  stored  for  the  further  solution. 


Having  found  W(^y)  for  the  rectangular  window  at  the  points  the  G  '(4)  of  equation 
(4.1 1)  may  now  be  interpolated  for  a  given  value  of  ^  =  /(^p^2)  within  the  window  as 


G-'(q)  =  C,{q„q,)  Wfe) 

where 


(5.1.11) 
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(5.1.12) 


CoiMl)  =  QA  ^1^2 

^  m  ml  ^IqI) 

now  p(q)  (equation  (4.5))  at  the  value  of  ^  =  fiq^A)  within  the  window  is  then  approxi¬ 


mated  as 


m  =  0-'{g)  H(q)  (5.1.13) 

and  so  equation  (4.8)  can  be  evaluated  at  q  in  order  to  find  the  reduced  system  operator 


Lo(^)  =  L„(^)+L,3(^)  m 

from  which  the  reduced  system  equation  can  be  solved  for  the  response  u„ 


(5.1.14) 


Lo  =  4 


(5.1.15) 


It  should  be  noted  that  the  matrix  )3( q)  is  a  COMPLEX  “Slave”  DOF  by  “Master”  DOF 
matrix  and  MUST  be  stored  if  the  “Slave”  response  is  to  be  reconstructed  as 


Us  =  ^(?)  u„ 


(5.1.16) 


5.2  Hermitian  Interpolation  of  ID  Linear  “Window” 

5.2.1  End  Point  Hermitian  Interpolation 

If  there  is  a  desire  to  confine  the  analysis  to  a  frequency  range  (line)  which  is  either  along 
the  “Real”  or  “Imaginary”  axes,  the  above  equations  may  be  greatly  simplified.  The  criteria  are 
for  a  line  along  the  “Real”  axis, 

G),  0  ,  Glj  =  0  ,  ^  0  ,  =0 

and  for  a  line  along  the  “Imaginary”  axis, 

G),  =  0  ,  G)2  ^  0  ,  -  Q  ,  0 
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Here  the  interpolation  points  of  a  frequency  range  or  line,  aj{a>i,a)2,w^,h^)  =  0)q,0)^ 
are  defined  as  the  line  end  points  ;  cDq  =  co^  +  i  co^  is  the  origin  or  beginning  of  the  line  and 

eitherco^  =  (to,  +  w^)  +  i  O)^  is  the  end  point  of  the  line  for  the  case  where 

^  ^  ’  K  =  the  “real”  axis  or  (0,  =  to,  +  i  (to^  +  h^)  for  the  case  where 

Wo,  =  0  ■,  ^  0,  i.e.  along  the  “imaginary”  axis. 


The  G"‘(<o)  is  now  simply  be  approximated  at  any  point  along  the  line, 
a  =  /(®„tb2)  ,as 


G-\d>)  =  l^(d))  +  to  W^iw)  +  d)^  W,{m)  +  d>'  W4(©)  (5.2.1) 

The  complex  constructor  matrices  Wj{m)  are  determined  similarly  to  the  2D  case,  for  the  line 
end  points,  o5j ,  as 


r  1 

dG  ’(<Oo) 

W2 

W3 

=  B-‘(tOo,to,)  0 

dco 

G-'K) 

3G-'(m,) 

y  do)  ) 

(5.2.2) 


and  the  real  interpolation  coefficient  matrix  B  *(£Oo,toJ  is  the  matrix  inverse  of 


5(tOo,toJ  =  B{co^,(02,w^,K)  = 


^1 

tOo 

0)1  ^ 

0 

1 

2ft)o 

3t0o 

1 

.0 

1 

20), 

30)1 

(5.2.3) 


As  in  the  2D  “window”  case  above,  the  Hermitian  interpolations  of  the  first  derivative  of  G  \q) 
with  respect  to  ^  are  as  above  in  equation  (5.1.1) 
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(5.2.4) 


im  G-(e) 

dqj  dqj 


ag~'(g;) 


where  first  derivative  of - may  be  obtained  from  equation  (5.1.4)  via  the  chain  rule  as 

dot): 


dG{q)  Mj  _  (dL^^jq)  _  \ 

dqj  dcOj  dq  dq  dq 


dqj  d(Dj 


(5.2.5) 


As  before  from  the  definition  ^  +i^2i.e.  ^  =  im,  the  respective  derivatives  are 


=  i  and  4^  =  -1  (5.2.6) 

d(0^  da^ 

The  cases  for  a  line  parallel  the  “Real”  axis,  (O^  ^  0  ,  Oj  ^  0  ,  =0 

or  for  a  line  para/ZeZ  the  “Imaginary”  axis,®,  ^0,0)2  ^  0  ,  =  0  ,  ?£  0  ,  the  fre¬ 

quency  can  be  defined  as  ®  =  1®|  e'  "  +  a  where  |®1  is  a  frequency  magnitude,  6  is  a  real 

fixed  value  of  either  zero  or  —  and  a  is  a  complex  origin  offset.  If  0)  parallel  to  either  the  real  or 

2 

imaginary  axis  then  a  ^  0.  So  for  (O  parallel  to  the  real  axis  for  0  =  0  and  a  —  0  +  ia2, 

7Z 

i.e.  (Oj  is  offset  from  the  real  axis  by  ia2orfor  ^  ~  ^  ^  offset 

from  the  imaginary  axis  by  a,.  Then  the  values  for  ®o  and  (O,  represent  magnitudes  as  defined 
as 


K1  = 

where  *  represents  the  complex  conjugate. 


(5.2.7) 


Note  for  the  cases  where  0)  is  on  an  axis  a  —  0 ;  if  0  —  0  ;  i.e.  0)  along  to  the  real  axis 
or  if  0  =  —  ;  i.e.  ©  along  to  the  imaginary  axis,  the  magnitude  is  simply 
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where  *  represents  the  complex  conjugate. 


5.2.2  Generalized  ID  Hermitian  Interpolation 

The  generalized  ID  Hermitian  Interpolation  allows  more  than  two  interpolation  (inver¬ 
sion)  points  and  can  utilize  higher  order  Hermitian  Interpolations.  However,  here,  we  shall  re¬ 
strict  the  analysis  to  1®'  order  Hermitian  Interpolation  with  up  to  4  interpolation  points.  Also  the 
equations  shall  be  derived  in  terms  of  the  generalized  variable  q  =  ico. 


The  G  '(4)  may  be  approximated  at  any  point  on  a  line,  q  =  f{q^,q2)  .  as 


G"(«)  =  W,(g)  +  «  W,(q)  +  ^  W,(q)  +  -  +  V/^^,„,(q)  (5.2.9) 


It  should  be  noted  that  for  simplicity,  interpolation  will  be  restricted  to  a  line  either  on  the  “Real” 
axis  or  a  line  on  the  “Imaginary”  axis.  So  from  the  definition, 

q  =  i(o=  i[(0^+i(02)  then  q  =  qi  +  iqi  =  -<^2  ■ 

Therefore,  for  real  frequencies,  rOi  ,  on  the  “Real”  frequency  axis 

m,  56  0  ,  <<>2  =  0  »  ^  0  ,  =0  and  so 

qi  =  0  ,  q^  ^  0  ,  =  0  ,  ^  0 

and  for  imaginary  frequencies,  0)2  .  the  “Imaginary”  frequency  axis 

=  0  ,  ©2  56  0  ,  =  0  ,  0  and  so 

qi  ^  0  ,  ^2  =  0  ,  ^  0  ,  =0 


(See  Appendix  B) 


For  r‘  order  Hermitian  Interpolation,  there  will  be  2N ,  '\.&.{Order  +  \)N ,  W  matrices 
needed  for  the  approximation  of  G~\q)  in  equation  (5.2.6).  These  matrices  may  be  calculated  as 


IT, 


=  B-\q)  o 


G-%)  ] 

dq 

G  *(^2) 

9G~*(^3) 

dq 


G-'fc) 

^G  ^(q^) 


(5.2.10) 


NOTE:  r‘  order  Hermitian  Interpolation  requires  only  the  first  derivative  of  G  \q)  with  respect 
to  q.  Higher  order  Hermitian  Interpolation  would  require  higher  order  derivatives. 


The  interpolation  coefficient  matrices,  B(q),  for  different  number  of  interpolation  points, 
N  =  2 , 3  ,  or  4,  and  for  f  ‘  order  Hermitian  Interpolations  may  be  written  explicitly  as 


For  N  =  2 


For  N  =  3 


B{q) 


^1 

Qi 

qx 

qx  ^ 

0 

1 

2qx 

^qx 

1 

% 

ql 

ql 

1 

'^qi 

nl) 

^1 

qx 

qx 

ql 

ql 

ql  ^ 

0 

1 

2% 

^ql 

^ql 

5^; 

1 

qi 

ql 

ql 

ql 

ql 

0 

1 

'^qi 

^ql 

4ql 

5ql 

1 

^3 

ql 

ql 

ql 

ql 

.0 

1 

2^3 

3q^ 

5qlj 

(5.2.11a) 


(5.2.11b) 
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For  N  =  4 


m 


^1 

Qt 

^1’  ^ 

0 

1 

2^1 

3q^ 

5q: 

et, 

1 

^2 

^2 

«2 

^2 

^2 

^2 

0 

1 

2^2 

'It 

1 

% 

^3 

^3 

0 

1 

2^3 

3ql 

5^3“ 

It 

1 

^4 

^4 

^4 

^4 

^4 

^4 

t 

.0 

1 

2^4 

4^4 

5^; 

6^4' 

'^tj 

(5.2.11c) 


As  above  in  the  2D  “window”  case  above,  the  Hermitian  interpolations  of  the  first  deriva¬ 
tive  of  with  respect  to  ^  are  as  above  in  equation  (5.1.1) 


dqj 


-G-(«)  ^  G-'W 

dqj 


(5.2.12) 


The  first  derivative  of  may  be  obtained  from  equation  (5.1.4)  and  from  the  definition 

dq 

of  a  =  <7,  +  ia,  with  the  derivatives  ^  ^ 

^  ‘  dqj  dq^  dq^ 

dG{q)  _  dL^{q)  _  ,  ^  dLi^jq)  _  da{q)  ^ 

^(£7..  da  dq  dq  )  dq: 


K 


(5.2.13) 


6  Lagrangian  Interpolation  Method 

6.1  2D  Interpolation  using  Rectangular  Lagrangian  Shape  Functions 

Isoparametric  2D  quadratic  shape  functions  of  a  rectangular  “element”  can  be  used  to 
interpolate  any  point  within  the  element.  A  rectangular  FEM  element  with  4  to  9  nodes  having 
the  following  node  numbering  is  used  to  represent  the  frequency  window  to  be  interpolated. 
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4 


7 


3 


The  2  dimensional  shape  functions,  hi{r,s),  in  isoparametric  coordinates  {r,s)  for  any 
nodal  configuration  are  defined  in  the  table  below. 

Table  6.1 

Include  only  if  Node  i  is  Defined 


i  =  5 

i  =  6 

i  =  7 

i  =  8 

i  =  9 

h, 

0.25  (1-r)  (1-5) 

-0.5  hs 

-0.5  hi 

-0.25  hg 

hz 

0.25  (1+r)  (1-5) 

-0.5  hs 

-0.5  h, 

-0.25  hg 

h. 

0.25  (1+r)  (1+5) 

-0.5  h. 

-0.5  hg 

-0.25  hg 

0.25  (1-r)  (1+5) 

-0.5  hj 

-0.5  hi 

-0.25  hg 

^5 

0.5  (l-A^)(l-5) 

-0.5  hg 

^6 

0.5  (l+r)(l-/) 

-0.5  hg 

h  j 

0.5  (l-r^)(l+5) 

-0.5  hg 

hi 

0.5  (l-r)(l-5') 

-0.5  hg 

hg 

(l-r^)(l-/) 

From  this  table,  the  2D  shape  functions,  h^{r,s),  for  the  represented  “element”  above,  that  is  the 
nine  nodal  points,  can  be  written  explicitly  as 
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h^{r,s)  =  (l-r^)  (l-/) 

h^{r,s)  =  0.5  (1-r)  -  0.5  h^{r,s) 

h,{r,s)  =  0.5  (l  -r^)  (l  -5)  -  0.5  h^{r,s) 

h^{r,s)  =  0.5  (1+r)  (l -  0.5  h^{r,s) 

h,{r,s)  =  0.5  (l-r")  (l -5)  -  0.5  h^{r,s)  (6.1.1) 

h^{r,s)  =  0.25  (1-r)  (I+5)  -  0.5  hj{r,s)  -  0.5  h^{r,s)  -  0.25  h^{r,s) 

hj{r,s)  =  0.25  (1+r)  (1 +5)  -  0.5  ^^(r,^)  -  0.5  hj{r,s)  -  0.25  hg{r,s) 

hi{r,s)  =  0.25  (1+r)  (I-5)  -  0.5  h^ir,s)  -  0.5  \{r,s)  -  0.25  h^{r,s) 

hi{r,s)  =  0.25  (1-r)  (1-s)  -  0.5  h^{r,s)  -  0.5  h^{r,s)  -  0.25  h^{r,s) 

Defining  ^1,^2  ^  the  origin  of  the  rectangular  window  to  be  interpolated,  as  the  win¬ 
dow  width  and  as  the  window  height,  the  nodal  point  coordinate  values  are  then  defined  as  in 
the  figure  below. 


w„ 


^l+'y’?2+\  <ll+^q^<l2+K 


h 

<lv^2+-r 


^1 


«1+>V,,92 


The  matrix  G{q)  as  in  equation  (4.6a)  and  (4.7a) 

G(?)  =  L22(9)-a(^)  L,2(^) 

=  (9^M22  +  4D22  +  ^22)~^(^)  +  ^D,2  +  8,2) 


(6.1.2) 
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must  be  evaluated  at  each  of  the  nine  nodal  points,  (or  fewer  if  less  interpolation  points  are 
used)  so  the  inverse  matrix  can  then  be  obtained  at  each  of  the  interpolation  points  and 

saved  for  further  use.  It  should  be  noted  that  there  will  be  a  maximum  of  9  COMPLEX  “Slave” 
DOF  by  “Slave”  DOF  matrices  of  the  respective  at  the  nine  nodal  points  stored  for  the 

further  solution. 


For  a  given  value  of  ^  =  /(4,42)  within  the  window  i.e.  <4  and 

q^  <  ^2  -  ^2  ^9’  the  corresponding  isoparametric  coordinates  {r,s)  can  be  found  from  the 

following  relationships 


-2  (0.5  -  gi) 

% 


(6.1.3a) 


^  ^  ^1 _ M  (6.1.3b) 

K 

The  value  of  G~\q)  can  be  interpolated  by  the  shape  functions  evaluated  at 
q  =  from  the  isoparametric  coordinates  (r,^)  and  the  matrices  previously 

calculated  at  each  of  the  “element”  or  window  interpolation  points  by 


g-'(9)  =  lAAr.S)  am 

l*=l 

The  matrix  ^{q)  is  then  approximated  as 

P{q)  =  G-\q)H{q)  (6.1.5) 

This  is  now  use  to  determine  the  reduced  operator  Lo(^)  as  before 

Lo(^)  =  L,(^)  + 1,2(9)  m)  (6  16) 

=  +  9D,,  +  S„)  +  (9^M,2  +  9D,2  +  S12)  fi{q) 

The  reduced  system  equation  can  be  solved,  as  before,  for  the  response  u^ 


24 


and  “Slave”  response  can  be  reconstructed  by 


u,  =  P{q)  u„ 


(6.1.7) 


>3^  The  matrix  is  a  COMPLEX  “Slave”  DOF  by  “Master”  DOF  matrix  and  MUST  be 
stored  if  the  “Slave”  response  is  to  be  reconstructed. 


6.2  Lagrangian  Interpolation  of  ID  Linear  “Window” 

6.2.1  Fixed  Point  Lagrangian  Shape  Functions 

Isoparametric  ID  shape  functions  of  a  linear  “element”  can  be  used  to  interpolate  any 
point  along  a  selected  line.  A  ID  FEM  linear  element  with  2  to  4  nodes  having  the  following 
node  numbering  is  used  to  represent  a  line  either  along  or  parallel  to  the  Real  or  the  Imagi¬ 
nary”  frequency  axes. 


The  1  dimensional  shape  functions,  hi{r),  in  isoparametric  coordinate  r  for  any  nodal  configu¬ 
ration  are  defined  in  the  table  below. 


Table  6.2 


Include  Only  if 
Node  3  is  Present 

Include  Only  if  Nodes 

3  and  4  are  Present 

K 

0.5  (  1  -  r ) 

-0.5 

+  .i(-9r^+r^  +  9r-l  ) 

/l2 

0.5  (1+r) 

-0.5  (l-r^) 

+  ^(9r^  +  ^-9r-l  ) 

h. 

(1-^) 

+  ±(27  r^  +  1  r^-21r-l  ) 

h4 

^(-27/^-9^  +  27r  +  9  ) 

NOTE:  For  the  3  node  shape  functions,  point  3  is  fixed  at  a  value  of  one-half.  For  the  4  node 
shape  functions,  point  3  is  fixed  at  one-third  and  point  4  is  fixed  at  two-thirds. 

From  Table  6.2,  the  ID  shape  functions,  /^,(r),  for  a  4  node  line  interpolation  ,  can  be 
written  explicitly  as 

Kir)  =  ±(-  27  -  9  +  27  r  -  9  ) 

h,{r)  =  (l-r^)  +  ^(  27  r'  +  7  r'  -  27  r  -  7  )  (6.2.1) 

h,{r)  =  0.5  (1+r)  -  0.5  (l-r^)  +  ±(  9  -f-  -  9  r  -  1  ) 

\{r)  =  0.5  (1-r)  -  0.5  (l-r")  +  ^(  -9  r'  +  r'  +  9  r  -  1  ) 

For  a  given  value  of  q  =  /{qiAi)  along  a  line  i.e.  <  ^,  +  %  if =  0  or 

q^  ^  ^2  -  ^2  K  ~  corresponding  isoparametric  coordinate  r  can  be  found  from 
the  relationship  in  equation  (6.1.3a)  as  either 


-2  (0.5  +  q,  -  q,) 


or 


-2  (0.5  K+  q2-  q2) 


ifh^  =  0 


if  =  0 


(6.2.2a) 


(6.2.2b) 
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The  value  of  can  be  interpolated  by  the  shape  functions  evaluated  at 

q  =  from  the  isoparametric  coordinate  f  and  the  matrices  previously  calcu¬ 

lated  at  each  of  the  line  interpolation  points  by 

i=l 

is  then  approximated  as  before  as  in  equation  (6.1.5) 


6.2.2  Variable  Point  Lagrangian  Shape  Functions 

If  it  is  desired  to  have  variable  locations  for  the  ID-shape  functions,  the  following  forms 
are  used. 


For  a  given  value  of  ^  =  /(4>42)  ^ilong  a  line  i.e.  q^  <  q  <q^  +  ,  if  /i,  =  0  or 

q^  <  q  ^  ^2  +  .  if  =  0.  the  term  Ax  is  defined  as 

Ax  =  ifh=0  (6.2.5a) 

% 

or 

Ax  =  ifw,  =  0  (6.2.5b) 
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Now  the  corresponding  3  point  quadratic  shape  functions  are  written  as 


Mq)  = 

,  X  +X,2 

X  Xj2 

Ax  + 

^(<?)  = 

J^12 

Ax  + 

X(x-Xi2) 

fh{q)  = 

X 

—  A  V  * 

^12  ”  ^12) 

ZjL/V 

- —  Ax^ 


X  X, 


12 

1 


X(x-Xi2) 

1 


Ax" 


(6.2.6) 


x,2(x-x,2) 


Ax" 


and  the  4  point  cubic  shape  functions  are 


j  ((^  Xi2)  +  (x  ^^13)  + (-^12  ^13))  ^  _  (  X.  Xj2  Xjj) 


Ax^  - 


X  X, 


12  -^13 


X  X12  ^13 


X  X,.  X 


Ax' 


12  -^13 


^12 


h^{q)  =  - 


_ ^  ^13 _ ^ _ _  Ax^ - = - 

^12(^“  ^12X^12  ”^13)  ^12X^12  ”^13)  ^12(^“^12X^~^13) 


X(X  -  X,2  )( -  X,3  )  X(X  -  JC12  )( -  J^,3  ) 


(-X-X13) 


Ax' 


(6.2.7) 


Ax' 


hM)  = 


X  X, 


12 


-Ax 


X  +  Xi 


12 


Xi3(x  X|3)(  Xj2+Xj3)  ^13)(-^12  *^13) 


Ax^  + 


Xi3(x-Xi3Xj^12-^I3) 


Ax' 


Here  x  is  defined  as  unity,  i.e.  x  =  1  and  then  the  terms  Xjj  and  X13  have  fractional  values  such 
that  0  <  X32  <  X33  <  1 . 

Equations  6.2.6  and  6.2.7  are  used  as  before  in  equation  6.2.3,  in  order  to  calculate  the 
value  of  G~\q)  at  q  =  tiq^Qi)  front  the  matrices  previously  calculated  at  each  of 

the  line  interpolation  points  by 

3or4 

G-(«)  =  Gr‘(?)  (6-2.8) 
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7  Results 


7.1  Introduction 

Initially,  the  FEM  models  were  generated  by  the  FWR  program  and  were  then  repro¬ 
duced  and  checked  against  identical  ABAQUS™  models.  ABAQUS™  was  then  used  to  develop 
more  complicated  FEM  models  for  use  in  the  FWR  method  (e.g.  spherical  shell  model).  That  is 
the  global  matrices  (mass  and  stiffness)  used  in  these  examples  were  generated  by  the  FWR 
program  for  the  simple  models  and  ABAQUS™  program  for  the  larger  more  complicated 
models.  The  matrices  were  then  reordered  and  partitioned  according  to  the  choices  of  “master” 
node  points. 


Calculations  were  preformed  using  both  the  Hermitian  and  Lagrangian  interpolation 
methods  in  order  to  obtain  a  reduced  solution  for  direct  comparison  to  solutions  of  the  full  DOF 
model.  The  full  DOF  solution  0^,,  is  found  from  the  equation 

Lu„  =  f  (7.U) 

And  the  reduced  solution  UReduced  is  calculated  from  the  reduced  DOF  equation 


Lo  u™ 


=  L 


And  the  reconstructed  “slave”  displacement  solution 


(7.1.2) 


(7.1.3) 


Hence,  the  total  reduced  solution  (“master”  and  reconstructed  “slave”  results)  is  then 


u 


Reduced 


^u  ^ 


(7.1.4) 


The  Hermitian  Interpolation  method  allows  2D  rectangular  window  interpolation  (4  p") 
and  ID  linear  interpolation  with  2  and  3or  or  variable  interior  point  location.  Note  that  the 
Hermitian  Interpolation  method  includes  the  function  interpolation  and  derivative  information, 
and  for  the  ID  cases,  with  3  or  4  points,  interior  window  information.  The  Lagrangian  Interpo- 
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lation  method  allows  2D  rectangular  window  interpolation  with  4,  5,  8  or  9  point  two  dimen¬ 
sional  shape  functions  and  ID  linear  interpolation  with  2 and  3  or  4  fixed  or  variable  interior 
point  location.  The  Lagrangian  Interpolation  method  for  2D  rectangular  interpolation  with  5,  8  or 
9  point  locations  and  ID  linear  3  or  4  points  have  interior  window  information. 

The  interpolation  method  for  the  2D  window  was  performed  on  several  of  the  following 
examples.  However,  since  the  frequency  sweeps  were  confined  along  the  real  frequency  axis,  the 
2D  analyses  gave  identical  results  when  compared  to  the  appropriate  ID  analyses.  The  2D 
Hermitian  interpolation  method  has  two  interpolation  (inversion)  points  on  the  real  axis,  hence 
its  results  were  the  same  as  the  ID  2'”  Hermitian  interpolation  method.  Similarly,  the  2D 
Lagrangian  method  which  has  three  interpolation  points  on  the  real  frequency  axis  was  identical 
to  the  ID  3*’'  Lagrangian  method. 

In  all  of  the  following  examples,  results  are  shown  for  the  ID  interpolation  method  con¬ 
fined  to  the  real  frequency  axis.  Also  the  analyses  for  all  examples  were  performed  using  the 
higher  point  interpolations,  4’’'  Hermitian  and  4'"  Lagrangian.  These  methods  gave  more  accurate 
results  for  comparison  to  the  full  DOF  solutions. 


7.2  Error  Analysis 

Error  analysis  was  performed  utilizing  the  relative  error  of  the  infinity  norms  of  the  ma¬ 
trices.  This  relative  error  is 


Relative  Error  = 


(M 

WMi 


(7.2.1) 


With  the  infinity  norm  of  A  is  defined  as  |1a|L  =  |  where  A  is  an  m  x  n  ma- 

;=i 

trix.  The  relative  error  as  defined  above  is  a  measure  of  the  decimal  digit  accuracy  of  A  as  com¬ 
pared  to  A. 
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From  the  definition  of  the  relative  error  equation  (7.2.1)  and  the  resultant  solutions  for  the 
full  DOF  problem  (via  7.1.1)  and  the  reduced  DOF  problem  (via  7. 1.2-7. 1.4),  the  relative  error  is 
then  written  as 


Relative  Error  = 


1  reduced 

^/b«) 

oo 

oo 

(7.2.2) 


Where  and  u^„  represent  the  entire  solution.  This  relative  error  represents  a  stringent 

global  error  estimate. 


The  magnitude  of  the  displacement  for  individual  DOFs  as  obtained  from  the  full  DOF 
solution  as  well  as  from  reduced  solutions  from  the  interpolation  methods  were  directly  com¬ 
pared.  Percent  error  analysis  was  performed  on  specific  DOFs  to  directly  compare  the  solutions. 
The  A%  or  the  percent  error  is  defined  as 

^reduced 

I  I 

Where  and  Uj^,i  here  represent  the  magnitude  Of  the  displacement  solution  of  a  specific 

individual  DOF. 


100 


(7.2.3) 


7.3  Eigenvalue  Analysis 

Eigenvalue  analysis  for  all  models  was  performed  with  ABAQUS™  and  LAPACK™ 
routines  for  inclusion  in  the  data  analyses.  The  eigenvalue  problem  can  be  written  in  two  ways, 
the  standard  eigenvalue  form  as 

A  X  =  Xx  (7.3.1) 

Or  the  generalized  eigenvalue  form  as 

Ax  =  A^Bx  (7.3.2) 
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The  frequency  dependent  mechanical  system  equation  is 

+  /<dD  +  S)m  =  /  (7.3.3) 

Where  M,  D  and  S  are  the  mass,  damping  and  stiffness  matrices  respectively.  The  term  u  is  the 
Fourier  transform,  F  =  e'“',  of  the  system  displacements  and  /  is  the  Fourier  transform  of  the 
system  loads. 

From  equation  (7.3.3),  if  the  mechanical  system  is  assumed  to  have  no  damping  term  pre¬ 
sent,  i.e.  (-£0^M  +  S)m  =  /,  then  the  generalized  eigenvalue  problem  is  written  via  equation 

(7.3.1)  as 

S  X  =  X^Mx  (7.3.4) 

Where  only  the  mass  and  stiffness  matrices  are  present. 

This  equation  may  be  solved  for  the  eigenvalues  directly  by  use  of  the  LAPACK™  rou¬ 
tines  LA_S YGV  or  LA_GEGV  for  symmetric  or  non-symmetric  matrices  respectively.  Alterna¬ 
tively  this  equation  may  be  solved  by  reformulation  via  the  Cholesky  decomposition  of  the  mass 
matrix,  M,  as 


M  =  UW  (7.3.5) 

The  decomposed  mass  matrix  is  then  used  to  transform  the  generalized  eigenvalue  problem 
(7.3.2)  into  the  standard  eigenvalue  problem  (7.3.1)  as 

Ax  =  vx  (7.3.6) 

Where  the  matrix  A  is 

A  =  (u^y'su  (7.3.7) 

And  the  eigenvalues  are  defined  as  v  =  A  ^  =  (2;r<o)^ .  The  eigenvalues  may  now  be  solved  for 
by  use  of  the  LAPACK™  routines  LA_POTRF,  LA_SYGST  and  LA_SYEV. 
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The  generalization  of  the  standard  eigenvalue  problem  to  nonlinear  eigenvalues  is 
(AAVBA  +  C)  x=  0  ('7-3.8) 


Which  can  be  seen  to  be  a  quadratic  equation  in  X.  This  nonlinear  problem  can  be  transformed  to 
a  2n  X  2n  “standard”  eigenvalue  problem  by 


"  0 

■  1 

fx^ 

=  A 

fx^ 

^-A"*  C 

T 

< 

1 

^yj 

Where  I  is  the  Identity  matrix,  the  matrix  A  is  defined  as 

f  0  I  ^ 

A  = 

l^-A-'C  -A-‘bJ 


And 


X  = 


^yj 


So  the  2n  x  2n  “standard”  eigenvalue  problem  now  becomes 


A 

A  3c  = 


vjc 


(7.3.9) 


(7.3.10) 


(7.3.11) 


(7.3.12) 


Applying  this  reformulation  to  the  mechanical  system  equation  (7.3.3)  when  damping  is 
present,  the  matrix  A  for  the  standard  eigenvalue  problem  becomes 


f  0  I 

A 

A  = 

And  then  the  eigenvalues  are  now  defined  as 
V  =  iX  =  iTvio) 


(7.3.13) 


(7.3.14) 
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NOTE:  For  this  standard  formulation  it  is  assumed  that  it  is  possible  to  obtain  the  inverse  of  the 
mass  matrix,  i.e.  M”'  exists. 


When  the  mass  matrix  in  not  invertable,  i.e.  there  are  zero  values  along  the  diagonal  of 
the  matrix,  it  is  possible  to  place  small  values  for  the  zero  terms  (several  orders  of  magnitude 
below  the  lowest  value)  on  the  diagonal.  The  presence  of  these  small  non-zero  diagonal  terms 
should  now  allow  the  inversion  of  the  mass  matrix  directly,  without  significantly  effecting  the 
accuracy. 

Alternatively  in  the  event  that  the  mass  matrix  is  not  invertable,  the  generalized  nonlinear 
eigenvalue  problem  (7.3.8),  i.e.  when  damping  is  present,  can  be  transformed  to  a  2n  x  2n  “gen¬ 
eralized”  eigenvalue  problem  by 

0  I 

=  X 

-C  1,0 

Or 

0  I  0 

=  X 

-S  — dJ  V^O  M  j  yyj 

The  generalized  eigenvalue  problem  is  then  becomes 

A  J  =  vBJ  (7.3.17) 

f  0  I  ^  n  0^ 

Where  the  terms  are  now  A  =  ,B  s  ,andv  =  iX  =  Ittico. 

I-S  -dJ  {o  m) 

NOTE:  For  both  the  standard  and  generalized  eigenvalue  problems  defined  above  for  the  nonlin¬ 
ear  case  2n  eigenvalues  are  obtained  which  will  result  in  DET  A  —  AI  =0  for  the  standard 

eigenvalue  problem  and  DET|a  -Ab|  =  0  for  the  generalized  eigenvalue  problem.  These 

eigenvalues  also  will  satisfy  DET|  (A  ^A  +  AB  +  C) |  =  0 . 


(7.3.15) 

(7.3.16) 
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Both  of  the  above  formulations  for  the  nonlinear  eigenvalue  problem  in  the  case  where 
damping  is  present  were  tested  via  the  LAPACK™  routines  LA_GEEV  and  LA_GEGV.  Both 
methods  resulted  in  identical  eigenvalues.  These  values  were  checked  for  consistency  by  apply¬ 
ing  the  appropriate  matrix  determinant  calculation. 


7.4  Program  Description 

All  computational  codes  were  developed  in-house  at  NRL  using  standard  FORTRAN  90 
(F90).  The  codes  were  initially  written  and  tested  on  an  Apple  Macintosh™  Power  PC  computer 
using  Apple  Macintosh  Programmers  Workshop™  (MPW)  and  Absoft  ProFortranT^^  version  6.0. 
These  codes  have  been  successfully  ported  to  a  Windows  based  PC  platform.  Fortran  90  was 
chosen  because  of  its  many  programming  features.  The  dynamic  memory  allocation  of  arrays, 
access  to  subarrays  directly  by  indexing,  usage  of  pointers  as  well  as  many  implicit  matrix  rou¬ 
tines  are  feature  which  make  F90  very  attractive  for  developing  the  necessary  computational 
codes  for  this  methodology.  Additionally  F90  offers  the  option  of  developing  parallel  codes  for 
High  Performance  Computing  (HPC)  platforms.  In  F90,  codes  can  easily  be  developed  in 
modular  form  with  the  “modules”  (.mod  file)  acting  similarly  to  linked  static  libraries.  This  al¬ 
lows  the  development,  testing  and  re-use  of  routines  in  multiple  programs.  Standard  routines, 
such  as  LU  decomposition,  Cholesky  decomposition,  eigenvalue  solvers,  etc.,  were  implemented 
from  the  LAPACK™  F90  version  2.0.  LAPACK  was  incorporated  in  order  to  minimize  the  need 
to  develop  standard  matrix  computational  routines,  but  yet  allow  broad  cross-platform 
compatibility. 

Brief  descriptions  of  the  major  codes  for  the  FWR  method  and  a  listing  of  the  working 
directory  and  of  program/routine  names  follow  in  this  section. 
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Calculation  of  Global  and  Partitioned  Matrices 


FEM_GIobals 

This  is  an  in-house  developed  code,  which  reads  a  discretized  model  file  and  calculates 
the  global  mass  and  stiffness  matrices  for  the  FEM  model.  These  full  global  matrices  may  either 
be  generated  as  un-ordered,  i.e.  numbered  as  in  the  original  FEM  model,  or  as  re-ordered  ac¬ 
cording  to  a  “master”  node  list.  Note  that  for  the  re-ordered  option  that  both  the  full  global  matri¬ 
ces  and  the  partitioned  matrices  are  saved  in  packed  binary  format.  Also  a  zero  damping  matrix 
is  output.  Because  of  the  limited  element  types  able  to  be  programmed,  this  program  was  initially 
used  in  order  to  test  the  subsequent  programs  for  simple  models. 


FEM_Globals 


.MDL 


Globals 
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Massll  Massl2  Mass21  Mass22 

-K 


Stiffnessll  Stiffnessl2  Stiffness21  Stiffness22 


Dampingll  Dampingl2  Damping21  Damping22 


AbaqusReader 

This  program  reads  the  output  file  from  ABAQUS™,  which  contains  the  computed  stiff¬ 
ness  and  mass  matrices  stored  in  lower  triangular  form.  The  program  generates  the  full  symmetry 
global  stiffness  and  mass  matrices  and  outputs  them  in  packed  binary  format.  Note  that  these 
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matrices  are  in  the  un-ordered  form,  i.e.  numbered  according  to  the  FEM  model.  A  zero  damping 
matrix  is  also  output. 


Stiffness  Mass  Damping 


Reorder_Globals 

This  program  applies  a  permutation  operation  on  the  full  global  mass,  stiffness,  and 
damping  matrices  based  on  a  “master”  node  list.  The  re-ordered  full  global  matrices  and  the  par¬ 
titioned  matrices  are  then  output  in  packed  binary  format. 


Massll  Massl2  Mass21  Mass22 


Stiffnessll  Stiffnessl2  Stiffness21  Stiffness22 


Dampingll  Dampingl2  Damping21  Damping22 
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RayleighDamping 

The  Rayleigh  damping  matrix  is  calculated  from  the  global  mass  and  stiffness  matrices 
by 


D  s  aM  +  iSS  (^Al) 

either  from  input  Rayleigh  damping  coefficients  Of  and  or  from  two  supplied  resonant 
frequencies,  O),  and  (o^,  and  their  respective  percentages  of  critical  damping,  ^and  4.  The 
Rayleigh  damping  coefficients  are  found  as 


_  2(0)24  ^i^\)  (7.4.2) 

col -col 

and 

a  =  2o)i^,  -  j8o)f  =  0),(2^,  -  fico^).  (7.4.3) 


Rayleigh  damping  may  either  be  calculated  for  unordered  or  re-ordered  matrices.  If  the 
calculation  is  for  the  unordered  matrices,  the  above  program,  “Reorder_Globals”  must  be  used  to 
obtain  the  re-ordered  full  global  matrices  and  the  partitioned  matrices  for  use  in  the  FWR  method 
programs. 


RayleighDamping 


Stiffness  Mass 


Damping 


OR 
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RayleighDamping 


Globals  I - K 


Stiffness  Mass 


Globals 
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Damping 
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Dampingll  Dampingl2  Damping21  Damping22 
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Calculation  of  Frequency  Window  Reduction  Method 


FWR_Constructors 

From  the  input  frequency  window  location  and  size,  the  interpolation  parameters  and  the 
saved  partitioned  global  matrices,  the  interpolation  constructor  matrices  may  be  calculated  for 
both  Lagrangian  and  Hermitian  interpolation  methods.  For  the  Lagrangian  interpolation  method 
(Section  6)  at  each  interpolation  point,  2  <  N  <  4,  the  inverse  of 

G{q)  =  t^J^q)-a{q)  L,2(^)  (7.4.4.) 

is  calculated  and  saved  in  packed  binary  format. 


The  interpolation  constructor  matrices,  for  the  Hermitian  interpolation  method  for 

a  ID  linear  interpolation  with  2  <  N  <  4  (Section  5)  are 


G-\q.) 

*(^i) 

dq 

dq 

G-\q^) 

dq 
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dG~\q^,q^) 

dq^ 

d^G-\q„q,) 

dq^dq^ 


«'(«/)  =  8"te)  “ 


G'\qi+w,q^-\-h) 

dG~\q^+w,q2-^h) 

dq^ 

dG~\qi+w,q2+h) 

dqi 

+w,g2  +^) 

dqidqt  ) 


(7.4.6) 


for  the  2D  rectangular  window  Hermitian  interpolation.  The  w[qj)  matrices  are  saved  a  packed 
binary  files  for  the  next  program  calculation. 


.RORD 


Globals 


FWR 


Constructors 


Results  Solution  Temporary  User 


t 
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FWR  LOs 


This  program  calculates  the  full  DOF  system  and/or  the  reduced  DOF  system  matrix  for  a 
given  interpolated  frequency  window  and  a  given  frequency  sweep  within  the  window.  Note  that 
the  appropriate  constructor  matrices  must  have  been  previously  calculated  by 
“FWR_Constructors”. 

Having  found  the  Hermitian  constructor  matrices,  ,  for  a  2D  rectangular  window 

at  the  points  qj,  may  be  interpolated  for  a  given  value  of  the  frequency  in  generalized 

coordinates^  =  within  the  window  as 

C-«)  =  W'fe)  (7-4.7) 

where 


=  (i  4  Qx  qi  Wi  qUi 

^  (7.4.8) 

Ql  ^1^2  ^1^2  ^2  ^1^2  ^1^2  ^1^2) 

Alternatively,  the  G~\q)  may  be  approximated  at  any  point  on  a  ID  line,  q  =  /(^i,^2) 

G-'(«)  =  W5(?)  +  q  WAq)  +  e  IV,®  +  -  +  s'"’"''"""!-'  (7.4.9) 

where  the  order  of  the  Hermitian  is  restricted  to  first  order  and  2  <  N  <  4. 


The  Lagrangian  methods,  having  calculated  the  appropriate,  G  *(^) ,  the  2D  interpolation 


IS 


G-'(9)  =  G-'(q) 

i=X 

and  the  ID  linear  interpolation  is 

G-\^  =  G-\q). 


(7.4.10) 


(7.4.11) 


i=l 
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The  number  of  interpolation  points  is  an  input  parameter  and  the  parametric  coordinates  are 
found  from  the  frequency  in  generalized  coordinates. 

The  having  calculated  G~\q),  the  matrix  f5{q)  is  given  as 

l}(q)  =  G-'{q)H{q)  (7.4.12) 

for  both  interpolation  methods  where  s  a(q)t,,(q)-  Lj|(j).  Finally  the  reduced  DOF 

system  constitutive  matrix  is 

Lo(4)  =  L„(4)+L„(^)  m  (7.4.13) 

Both  Lo(^)  and  i3(^)are  saved  in  one  packed  binary  format  file.  Also  the  full  DOF  system  ma¬ 
trix  may  be  calculated  and  saved  from, 

Lo  =  +  S  (7.4.14) 


FWR_LOs 


.FWR 

.RORD 


Constructor  Files 


FWR 


Constructors  Results  Solution  Temporary  User 


LO### 


t 


AND/OR 
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Solve_FWR 


With  the  data  calculated  above,  this  program  solves  the  full  and  reduced  DOF  system 
matrices  from,  Lj,  which  are  saved  as  a  function  of  frequency  across  an  interpolation  window. 

This  program  now  solves  for  a  given  load  and  boundary  conditions  the  equation  Lj  u  =  f  for  the 
full  DOF  system  and/or  the  equation  L,,  u^  =  f„  for  the  reduced  DOF  system.  Also  for  the  re¬ 
duced  system,  the  total  reconstructed  solution  is  calculated  from  u^  =  u^ .  The  solutions  as  a 
function  of  frequency,  individual  DOF  solutions  as  well  as  the  error  analysis  described  in  section 
7.2  are  all  outputs  of  this  program. 


LO  Files 
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FWR.EVs 


This  program  calculates  the  constrained  or  unconstrained  eigenvalues  of  the  full  DOF 
model  according  to  the  procedure  described  in  section  7.3  and  LAPACK  routines  [12].  Note  that 
the  program  will  calculate  the  eigenvalues  for  systems  with  and  without  damping. 


FWR_EVs  Output 


.SLVR 


Stiffness  Mass  Damping 
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Application  Directory  and  Working  Directories 


Globals 


Files  of  full  and  partitioned  Mass,  Stiffness  &  Damping  Matrices 


Unordered  piles  of  full  Unordered  Global  Matrices 


Models 


Solvers 


Contains  “NameList”  text  files  of  the  FEM  Model  Files 


Contains  “NameLisf  ’  text  files  of  the  FEM  Element  Types 

Elements 

Contains  “NameList”  text  files  of  the  Solver  Files 


Modules  Drivers 


Mains 


vlodule.  Driver  and  Main  Source  Text  Files 


Objs 


Contains  Compiled  Object  Files,  ‘‘.o”  and  “.mod’ 


FWR 


Storage  Directories  for  Frequency  Window  Reduction  Calculations 


ifilj  Ifllt  lili  1111 

Constructors  Results  Solution  Temporary  User 


LAPACK  Modules  and  Pre-Compiled  Libraries 


LAPACKMods  LAPACKLibs 


Directory  for  Storage  of  Frequency  Window  Reduction  Output  Files 


Run 
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Listing  of  Program  and  Routine  Names 


MAINS: 

PROGRAM  FEM_Globals 
PROGRAM  AbaqusReader 
PROGRAM  Reorder_Globals 
PROGRAM  RayleighDamping 
PROGRAM  FWR_Contructors 
PROGRAM  FWR_LOs 
PROGRAM  Solve_FWR 
PROGRAM  FWR_EVs 

DRIVERS: 

SUBROUTINE  Do_Global_Matrices 
SUBROUTINE  Do_Reorder_Global 
SUBROUTINE  Do_FWR_Contructors 
SUBROUTINE  Do_FWR_L 
SUBROUTINE  Do_Solve_FWR 
SUBROUTINE  Do_FWR_EV 

MODULES: 


MODULE  nrtype 
MODULE  FEM_BOUNDS 
MODULE  FEM_CONTROL 
MODULE  FEM_INC 
MODULE  FEM_INC_FIXED 


MODULE  standard_file_names 

FUNCTION  FILE_NAME(IO,ID) 


MODULE  sparseutil 

FUNCTION  CHECK_FILE(ID,IO_NAME) 

FUNCTION  CHECK_SIZES(M,N,MP,NP) 

FUNCTION  CHECK_MATMUL(M,N, message) 

SUBROUTINE  CHECK_SPARSE(IO_NAME,len,m,n) 
SUBROUTINE  CHECK_THRESHHOLD(A,threshhold, OPTION) 
FUNCTION  arth(first,increment,n) 

SUBROUTINE  GET_SPARSE(IO_NAME,spa) 

SUBROUTINE  PUT_SPARSE(IO_NAME,spa) 

SUBROUTINE  sparsein(a,threshhold,spa) 

SUBROUTINE  sparseout(spa,threshhold,a) 

FUNCTION  matrixout(spa,m,n) 
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SUBROUTINE  PUT_MATRIX(IO,ID,MP,NP, matrix) 
FUNCTION  GET_MATRIX(IO,ID,MP,NP) 


MODULE  standard_interface 

SUBROUTINE  GET_DATE_TIME(ID) 

FUNCTION  CHECK_GLOBALS(N,OPTION) 

SUBROUTINE  CHECK_GLOBALS_MESSAGE(TESTFLAG,OPTION) 
SUBROUTINE  PRINT_GLOBALS_MESSAGE(IO, ID) 

FUNCTION  OVERWRITE_MESSAGE(IOPT) 

SUBROUTINE 

CHECK_CONSTRUCTORS(w,h,HNUM,SFNUM,HExist,SFExist) 

SUBROUTINE 

CHECK_CONSTRUCTORS_MESSAGE(TESTFLAGl,TESTFLAG2) 
SUBROUTINE  CONSTRUCTORS_OPTION(ISTART, BEND) 
SUBROUTINE 

FWR_L_OPTION(TESTFLAGl,TESTFLAG2,ISTART,IEND) 
SUBROUTINE  GET_FEM_FILE(ID,IO_NAME) 

SUBROUTINE 

CHECK_THRESHHOLD_OPTION_d(MP,NP,matrix,MATRIX_NAME) 

SUBROUTINE 

CHECK_THRESHHOLD_OPTION_z(MP,NP, matrix, MATRIX_NAME) 
SUBROUTINE  WRITE_MATRIX(X,MP,NP,M,N) 

SUBROUTINE 

PRINT_MATRIX_OPTION_d(MP,NP, matrix, MATRIX_NAME) 
SUBROUTINE 

PRINT_MATRIX_OPTION_z(MP,NP,matrix,MATRIX_NAME) 
FUNCTION  TRACE_EXECUTION(IOPT) 

MODULE  MTRXUTIL 

FUNCTION  FACTORIAL(N) 

FUNCTION  LA_LU_DETERMINANT(n,a) 

FUNCTION  LA_LU_INVERSE(n,a) 

FUNCTION  LA_LU_SOLVE(n,a,b,ID) 

FUNCTION  IS_MATRIX_SQUARE(A) 

FUNCTION  SYMMETRY_CHECK(A,threshhold) 

SUBROUTINE  SORTER(N,A,B,C,ID) 

FUNCTION  MATRIX_SYMMETRY(N,A) 

FUNCTION  IDENTITY_MATRIX(N,X) 

FUNCTION  GET_MATRIX_DIAGONAL(A) 

SUBROUTINE  PUT_MATRIX_DIAGONAL(A,X,THRESHHOLD) 
FUNCTION  ROW_COLUMN_ZERO_CHECK(N,A,ID) 

SUBROUTINE  QUERY_MATRIX(A,OPTION) 

FUNCTION  MATRIX_POWER(K,N,X) 

FUNCTION  A_MATRIX_ms(m,n,mass,stiffness) 

FUNCTION  A_MATRIX_msd(m,n,mass,stiffness,damping) 

SUBROUTINE  EIGEN_VALUES_s(m,stiffness,IsSSym) 
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SUBROUTINE  EIGEN_VALUES_ms(in,mass,stiffness,IsMSym,IsSSym) 
SUBROUTINE  EIGEN_VALUES_msd(m, mass, stiffness, damping,  & 
IsMSym,IsSSym,IsDSym) 

SUBROUTINE  WRITE_EIGEN_VALUES(ID,K,m,WR,WI,BETA,IO) 

MODULE  FWRUTIL 

FUNCTION  Complex_Magnitude(x,a) 

SUBROUTINE  CHECK_Q_BOUNDS(IO,wnd,pnt,InWnd) 

FUNCTION  GET_POINT_2D(ID,xO,w,h) 

FUNCTION  GET_POINT_lD(ID,xO,w,h,x  1 2,x  1 3) 

FUNCTION  dxi(ID,x,xO,w,h) 

FUNCTION  L(ID,MP,NP,x,derivative,option) 

FUNCTION  G(ID,MP,NP,x,derivative,option) 

FUNCTION  H(ID,MP,NP,x,option) 

FUNCTION  C(N,NUM,xl,x2) 

FUNCTION  B(N,NUM,xO,w,h,xl2,xl3) 

FUNCTION  C0(N,NUM,xl,x2) 

SUBROUTINE 

U_MATRICES(ID,MP,NP,NUM,xO,w,h,xl2,xl3,length,m,n,option) 
SUBROUTINE  W_M  ATRICES(ID,MP,NP,NUM,xO,w,h,x  1 2,x  1 3, option) 
SUBROUTINE  G_INVERESES(ID,MP,NP,NUM,xO,w,h,xl2,xl3,option) 
FUNCTION  APPROX_G_INVERSE(ID,IO,NUM,MP,NP,x,xO,w,h,xl2,xl3) 
FUNCTION  INTERPOLATION_FUNCTIONS_2D(M,MP,x,xO,w,h) 

FUNCTION 

INTERPOLATION_FUNCTIONS_lD(M,MP,x,xO,w,h,xl2,xl3) 

FUNCTION 

BETA_MATRIX(ID,IO,NUM,MP,NP,x,xO,wx,hx,xl2,xl3,option) 
FUNCTION  LO_MATRIX(ID,IO,NUM,MP,NP,x,xO,w,h,xl2,xl3 
,BETAMATRIX,option) 

FUNCTION  REDUCED_FWR_MATRICES(IO,MP,NP,BETAMATRIX) 

MODULE  USERUTIL 

FUNCTION  ALFA(ID,MP,NP) 

FUNCTION  ALPHA(ID,MP,NP,q) 

FUNCTION  dALPHA(ID,MP,NP,q) 

FUNCTION  ddALPHA(ID,MP,NP,q) 

FUNCTION  DIRICHLET_FUNC(T) 

FUNCTION  LOAD_FUNC(N,T) 

MODULE  GLOBALUTIL 

FUNCTION  ELASTIC_MATRDC(ID,NP,MP,IPRINT,PROPERTY) 
SUBROUTINE  INTEGRATION_PTS(IORDER,QUAD,RPT,SPT,TPT, 
RWT,SWT,TWT) 

FUNCTION  LOCAL_COORD(N,MP,NP,IPRINT,XYZ,CONNECT) 
SUBROUTINE  JACOBIAN(ID,M,N,IPRINT,XYZL,DSF,XJINV  ,DET) 
FUNCTION 
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B_MATRIX(ID,M,N,MP,IDCMTRX,IPRINT,XYZL,SF,DSF,XJINV) 
FUNCTION  H_MATRIX(ID,M,N,IDCMTRX,IPRINT,SF) 

SUBROUTINE  LOCAL_to_GLOBAL(NUMDOF,M,N,LP,NP,CONNECT, 
local,global,LOOKUP,OFFSET) 

SUBROUTINE  REORDER_GLOBAL(NDOF,N,NP, LOOKUP, 
global, reorder,NNDOF,ORGNLOFFSET,REORDOFFSET) 
SUBROUTINE  SHAPE_FUNCTIONS(ID,M,N,R,S,T,SF,DSF) 
SUBROUTINE  BAR_BEAM(M,N,MP,NP,R,S,T,SF,DSF) 

SUBROUTINE  RECTANGULAR(M,N,MP,NP,R,S,T,SF,DSF) 
SUBROUTINE  TRIANGULAR(M,N,MP,NP,R,S,T,SF,DSF) 
SUBROUTINE 

GLOBALS_OUTPUT(IO,NP,MP,MTRXTHSD,global,RENUMBER) 
SUBROUTINE  GLOBAL_SUMMARY(IO,ID) 

FUNCTION  OFF_SET_LIST(N,DOF_PER_NODE,LOOKUP) 
SUBROUTINE  MASTER_NODE_LIST(N,M,NNMF,NNFM) 

FUNCTION  QUAD_CHECK(IDELMT) 

FUNCTION  GET_SIDE_NUM(IDELMT,NUMNEL) 

SUBROUTINE  READ_ELEMENT_TYPE(ID,EL_INFO,STF_INT, 
MSS_INT,DMP_INT,FRC_INT) 

SUBROUTINE  PRINT_ELEMENT_TYPE(EL_INF0,STF_INT, 
MSS_INT,DMP_INT,FRC_INT) 

SUBROUTINE  SET_ELEMENT_SIDE(EL_INFO) 


MODULE  GLOBALS 

SUBROUTINE  GLOBAL_MATRICES 
SUBROUTINE  STIFFNESS(NO,NP,MP, global, extra,IDLIST) 
SUBROUTINE  MASS(NO,NP,MP,global,extra,IDLIST) 
SUBROUTINE  DAMPING(NO,NP,MP,global,extra,IDLIST) 

MODULE  MODELUTIL 

SUBROUTINE  ID_FEM_FILE(IO_NAME) 

SUBROUTINE  CLEAR_ALLOCATED_MATRICES(I) 
SUBROUTINE  FEMMODEL_FILE_IO(IO,IO_NAME) 
SUBROUTINE  FEMMODEL_SUMMARY(IO) 
SUBROUTINE  FWRMODEL_FILE_IO(IO,IO_NAME) 
SUBROUTINE  FWRMODEL_SUMMARY(IO) 
SUBROUTINE  CMRMODEL_FILE_IO(IO,IO_NAME) 
SUBROUTINE  CMRMODEL_SUMMARY(IO) 
SUBROUTINE  FEMSOLVER_FILE_IOaO,IO_NAME) 
SUBROUTINE  FEMSOLVER_SUMMARY(IO) 
SUBROUTINE  REORDER_FILE_IO(IO,IO_NAME) 
SUBROUTINE  REORDER_SUMMARY(IO) 

MODULE  SOLVERUTIL 

FUNCTION  BCMATRIX_d(MP,NP, matrix) 
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FUNCTION  BCMATRIX_z(MP,NP,matrix) 

FUNCTION  RHSFORCES(NP,T) 

FUNCTION  RHSBCS(NP,T, FORCES) 

SUBROUTINE  RUNGEKUTTA(MP,LP,NP,A1 3 1  ,AINV, 
NODE_NUM,NODE_DOF) 

SUBROUTINE  CENTRALDIFFERENCE(MP,LP,NP,A1 3 1  ,AINV, 
NODE_NUM,NODE_DOF) 

SUBROUTINE  NEWMARKSMETH0D(MP,LP,NP,A13 1  ,AINV, 
PRETERMS,NODE_NUM,NODE_DOF) 

SUBROUTINE  RECONSTRUCTION(MP,LP,NP,F,U,V) 
SUBROUTINE  RUN_SOLVERS(NODE_NUM,NODE_DOF) 
SUBROUTINE  SOLVER_PREP(MP,  A1 3 1  ,AINV, PRETERMS) 
SUBROUTINE  EXTRACT_TUV(NODE_NUM,NODE_DOF) 
SUBROUTINE  WRITE_ERROR_HEADER(I,ID,IO,NORM) 
SUBROUTINE  WRITE_RESULTS(ID, 10, LP,NDOF,LHS,norm_id) 
SUBROUTINE  WRITE_ERRORS(ID,LP,NDOF,ID2,norm_id) 
FUNCTION  CALCULATE_ERRORSa) 

SUBROUTINE  WRITE_ERROR_NORMS(ID,LP,ID2,norm_id) 
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8  Examples 

8.1  3-Layer  Steel-Polystyrene  ID  Model 

A  simple  one-dimensional  (ID),  3 -layer  laminate  was  considered  as  a  first  test  of  the 
FWR  method.  The  laminate  was  comprised  of  layers  of  steel  and  polystyrene  as  shown  in  Figure 

8.1.1  with  the  properties  listed  in  Table  8.1.1.  The  FEM  model  was  formed  from  26  ABAQUS™ 
2D  plain  stress  elements  consisting  of  46  node  points  representing  92  DOFs  (2  degrees  of  free¬ 
dom  per  node  point).  Fourteen  “master”  nodes  (28  DOFs)  were  selected  and  are  shown  in  the 
figure  as  symbols.  These  were  chosen  as  master  nodes  as  they  have  applied  loads  and  fixed 

boundary  conditions  in  the  FEM  model.  The  boundaries  between  element  groups  of  differing 
materials  and  element  grouping  midpoints  were  also  chosen  as  master  nodes.  This  represents  a 
70%  reduction  in  the  degrees  of  freedom,  i.e.  92  DOFs  for  the  full  problem  to  28  DOFs  for  the 
reduced  problem.  For  the  level  of  discretization,  only  frequencies  up  to  350  KHz  are  appropriate 
(assuming  8-10  elements  are  necessary  to  model  the  highest  frequency). 

Table  8.1.1  -  Properties  for  Steel-Polystyrene  Laminate 


Density 

Young’s  Modulus 

Poisson  Ratio 

Wave  Velocity 

Steel 

7900  Kg/m' 

256.7  GPa 

0.25 

5700  m/sec 

Polystyrene 

1050  Kg/m' 

9.387  GP  a 

0.25 

2990  m/sec 

The  largest  valid  frequency  window  of  0  to  350  KHz  was  examined  with  a  70%  reduction 
in  the  degrees  of  freedom.  Figure  8.1.2  shows  the  relative  error  for  a  frequency  sweep  over  the 
entire  frequency  window  of  0  to  350  KHz.  The  figure  compares  several  different  interpolation 
functions.  It  should  be  noted  that  at  the  end  points  and  the  two  interior  points  (approximately 
1 16.7  KHz  and  233.3  KHz)  the  error  is  very  small.  This  is  because  these  points  represent  points 
at  which  the  direct  inversion  for  the  solution  is  obtained;  hence  the  solution  is  nearly  exact.  It  can 
be  seen  that  generally  the  4^*  Lagrangian,  because  of  the  interior  point  information  is  better  than 
the  2P'  Hermitian,  even  though  the  Hermitian  has  derivative  information  at  the  interpolation  (end) 
points.  Clearly  the  4^'  Hermitian,  with  both  interior  point  function  and  derivative  information  has 
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the  lowest  relative  error  as  well  as  capturing  more  accurately  several  of  the  eigenvalues  of  the 
system.  Table  8.1.2.  The  eigenvalues  are  represented  on  the  plots  as  gray  vertical  lines.  Not  all 
the  eigenvalues  are  captured  because  of  the  frequency  sweep  step  size  of  1  KHz  is  not  fine 
enough  to  approach  the  singular  eigenvalue. 


Table  8.1.2  -  Calculated  Constrained  Eigenvalues 


Frequency  KHz 
0.909 

7.666 
19.522 
23.211 
61.726 

71.666 
119.722 
143.067 
167.134 
182.101 
256.836 
266.268 
300.511 
312.869 
318.306 
407.034 


The  actual  %  error  of  the  magnitude  of  the  resultant  displacement  for  the  “master”  DOF 
Ui,,  is  show  in  Figure  8.1.3.  The  4p'  Lagrangian  and  the  2'"  Hermitian  interpolation  functions  show 
very  large  errors  as  compared  to  the  full  DOF  solution.  However,  the  4*”  Hermitian  method  has 
generally  less  than  1%  error  over  the  entire  frequency  range.  The  comparison  of  the  4'" 
Lagrangian  and  the  4'*'  Hermitian  interpolation  to  the  full  DOF  solution  for  the  magnitude  of  the 
displacement  for  the  “master”  DOF  u,,(  is  show  in  Figure  8.1.4.  It  can  be  seen  that  that  both  pro¬ 
duce  reasonably  good  agreement.  The  %  error  of  the  magnitude  of  the  “reconstructed”  displace¬ 
ment  for  the  “slave”  DOF  Ujj^  shown  in  Figure  8.1.5,  as  expected,  is  identical  to  the  %  error  for 
“master”  DOFs  since  the  “slave”  DOFs  are  reconstructed  directly  from  the  “master”  DOF  solu¬ 
tion  according  to  equation  2.14. 

Using  the  frequency  window  range  of  0  to  350  KHz  for  the  interpolation  functions  and 
only  performing  a  partial  sweep  of  the  frequency  window,  it  is  possible  to  see  more  clearly  the 
effects  the  different  interpolation  functions  have  near  an  eigenvalue.  The  relative  error  for  an  in¬ 
terpolation  window  of  0  to  350  KHz  but  for  a  frequency  sweep  of  250  to  280  KHz  with  a  step 
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size  of  150  Hz  is  shown  in  Figure  8.1.6.  It  is  easily  seen  that  the  4^'  Lagrangian  interpolation  un¬ 
dershoots  the  eigenvalue  where  as  the  2*”  Hermitian  overshoots  the  eigenvalue  and  the  4'’* 
Hermitian  more  closely  matches  the  eigenvalue  response.  This  is  more  easily  demonstrated  in 
Figure  8.1.7  examining  the  magnitude  of  the  displacement  for  the  “master”  DOF  u,^which  shows 
how  closely  the  respective  solutions  approach  the  full  solution  near  this  eigenvalue.  These 
figures  indicate  how  accurate  or  smooth  the  interpolation  functions  are  within  a  given  interpola¬ 
tion  window. 

In  order  to  improve  the  accuracy  of  the  FWR  method,  the  frequency  window  interpola¬ 
tion  range  was  reduced.  This  reduction  in  the  window  range  allows  for  improvement  of  the  inter¬ 
polation  function  as  well  as  a  finer  frequency  sweep  step  size.  Also  with  a  reduced  frequency 
window  range,  there  are  fewer  eigenvalues  within  the  window,  which  improves  the  overall  accu¬ 
racy.  Figure  8.1.8  shows  the  relative  error  for  the  frequency  window  of  0  to  100  KHz  where  the 
sweep  step  size  was  500  Hz.  Dramatic  improvements  (several  orders  of  magnitude)  are  seen  in 
all  the  interpolation  methods.  This  dramatic  improvement  is  reflected  in  the  %  error  for  the  mag¬ 
nitude  of  the  displacement  for  the  “master”  DOF  shown  in  Figure  8.1.9.  Both  4^'  Lagrangian 
and  the  2”'  Hermitian  interpolation  methods  now  show  less  than  1%  error  and  the  4”*  Hermitian 
interpolation  method  is  less  that  0.00001  %  error.  The  plot  of  the  actual  the  magnitude  of  the  dis¬ 
placement  for  the  “master”  DOF  u,,^  in  Figure  8.1.10  shows  virtually  no  distinction  between  the 
interpolated  reduced  solutions  and  the  full  DOF  solution. 

The  further  reduction  of  the  interpolation  frequency  window  to  the  range  of  0  to  45  KHz 
with  a  frequency  sweep  step  size  of  150  Hz  is  shown  in  Figure  8.1.11.  Again  dramatic  improve¬ 
ments  (several  orders  of  magnitude)  are  seen  in  all  the  interpolation  methods  and  it  is  also  seen 
in  the  %  error  of  the  magnitude  of  the  displacement  for  the  “master”  DOF  u,^  shown  in  Figure 
8.1.12.  Here  the  Lagrangian  and  the  2**'  Hermitian  interpolation  function  are  less  than  0.1% 
error  with  the  d*"*  Hermitian  interpolation  method  is  less  than  0.0000001  %  error.  Clearly  for  this 
frequency  window  range  for  all  the  interpolation  functions,  the  magnitude  of  the  displacement 
for  the  “master”  DOF  u,,  is  numerically  identical  to  the  full  DOF  solution.  Figure  8.1.13. 
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For  the  interpolation  frequency  window  range  of  0  to  45  KHz,  the  number  of  DOFs  was 
dramatically  reduced  from  92  DOFs  to  8  DOFs  representing  only  the  loaded  and  constrained 
node  points.  This  represents  a  91%  DOF  reduction.  Comparing  the  relative  error  for  only  the  4^ 
Lagrangian  and  the  4'*'  Hermitian  interpolation  methods,  it  is  seen  in  Figure  8.1.14  that  despite 
such  a  drastic  reduction  in  the  number  of  degrees  of  freedom,  very  accurate  results  are  obtained 
for  this  model.  Even  with  a  91%  DOF  reduction  the  4^'  Lagrangian  interpolation  shows  approxi¬ 
mately  less  than  1%  error  and  the  4^'  Hermitian  interpolation  method  is  less  that  0.001  %  error  in 
Figure  8.1.15.  Also  the  actual  magnitudes  of  the  displacement  for  these  interpolation  methods  as 
compared  to  the  full  DOF  solution  are  nearly  identical  as  shown  in  Figure  8.1.16. 
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3-Layer  Steel-Polystyrene  ID  Figures 
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3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 


o 


(  jojja  )  OOl 


Figure  8.1.2-  Relative  Error  comparing  interpolation  functions  for  the  frequency  range  of  0  to  350  KHz 
and  a  frequency  step  size  of  1  KHz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Master  DOF  u^  ^  Displacement  Magnitude  %  Error 


Figure  8.1.3  -  %  Error  for  the  “master”  DOF  lulj^  (magnitude  of  the  displacement)  comparing  interpolation 
functions  for  the  frequency  range  of  0  to  350  KHz  and  a  frequency  step  size  of  1  KHz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Master  DOF  Displacement  Magnitude 
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Figure  8.1.4  -  The  “master”  DOF  |u|ix  (magnitude  of  the  displacement)  comparing  the  interpolation  functions, 
4ptLagrangian  and  Hermitian,  to  the  full  DOF  solution  for  the  frequency  range  of  0  to  350  KHz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Slave  DOF  u Displacement  Magnitude 


Figure  8.1.5  -  The  “slave”  DOF  [ujisx  (magnitude  of  the  displacement)  comparing  two  interpolation 
functions  to  the  full  DOF  solution  for  a  frequency  range  of  0  to  350  KHz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Interpolation  Window  0  to  350  KHz 
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Figure  8.1.6  -  Relative  Error  comparing  different  interpolation  functions  for  the  interpolation  window 
range  of  0  to  350  KHz  and  a  frequency  sweep  of  250  to  280  KHz  for  a  frequency  step  size  of  150  Hz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70  %  Reduction) 
Interpolation  Window  0  to  350  KHz 


Figure  8.1.7  -  The  “master”  DOF  |u|ix  (magnitude  of  the  displacement)  comparing  interpolation  functions 
for  the  window  of  0  to  350  KHz  and  a  frequency  sweep  of 250  to  280  KHz  at  a  step  size  of  1 50  Hz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 


Figure  8.1.8  -  Relative  Error  comparing  interpolation  functions  for  the  frequency  range  of  0  to  100  KHz 
and  a  frequency  step  size  of  500  Hz. 


Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Master  DOF  Ui„  Displacement  Magnitude  %  Error 


Figure  8.1.9  -  %  Error  for  the  “master”  DOF  |u|ix  comparing  interpolation  functions  for  a  frequency  range 
of  0  to  100  KHz  at  a  frequency  step  size  of  500  Hz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Master  DOF  u  Displacement  Magnitude 
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Figure  8.1.10  -  The  “master”  DOF  |u|ix  (magnitude  of  the  displacement)  comparing  the  4pt  Lagrangian  and 
the  4pt  Hermitian  interpolation  functions  to  the  full  DOF  solution  for  the  frequency  range  of  0  to  1 00  KHz. 


3  Layer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
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Figure  8.1.12  -  %  Error  for  the  “master”  DOF  luhx  comparing  interpolation  functions  for  a  frequency  range 
of  0  to  45  KHz  at  a  frequency  step  size  of  150  Hz.  NOTE:  70%  DOF  Reduction. 


yer  Laminate  :  92  DOF  to  28  DOF  (70%  Reduction) 
Master  DOF  u  Displacement  Magnitude 
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Figure  8.1.13  -  The  “master”  DOF  jujix  (magnitude  of  the  displacement)  comparing  two  interpolation 
functions  to  the  full  DOF  solution  for  a  frequency  range  of  0  to  45  KHz.  NOTE:  70%  DOF  Reduction. 


3  Layer  Laminate :  92  DOF  to  8  DOF  (91  %  Reduction) 


Figure  8.1.14  -  Relative  Error  comparing  the  4pt  Lagrangian  and  the  4pt  Hermitian  interpolation  methods 
for  the  frequency  range  of  0  to  45  KHz  and  a  frequency  step  size  of  150  Hz.  NOTE:  91%  DOF  Reduction. 


3  Layer  Laminate :  92  DOF  to  8  DOF  (91%  Reduction) 
Master  DOF  Ui„  Displacement  Magnitude 


Figvire  8.1.16  -  The  “master”  DOF  |u|ix  comparing  the  interpolation  methods  for  the  frequency  range  of  0 
to  45  KHz  and  a  frequency  step  size  of  150  Hz.  NOTE:  91%  DOF  Reduction. 


8.2  3-Layer  Steel-Polystyrene  2D  Model 

A  two-dimensional  (2D)  model  based  on  the  above  ID  model  of  the  3-layer  laminate  was 
considered  next.  The  FEM  model  was  formed  from  88  ABAQUS™  2D  plain  stress  elements 
consisting  of  1 15  node  points  representing  230  DOFs  (2  degrees  of  freedom  per  node  point) 
Figure  8.2.1.  Thirty-five  “master”  nodes  (70  DOFs)  were  chosen,  as  above,  and  are  shown  in  the 
figure  as  symbols.  This  represents  a  70%  reduetion  in  the  degrees  of  freedom,  i.e.  230  DOFs 

for  the  full  problem  to  70  DOFs. 

Initially  an  interpolation  frequency  window  of  0  to  25  KHz  was  examined,  however  it 
was  found  that  this  was  too  large  of  a  frequency  range  for  this  model.  A  smaller  interpolation 
window  of  0  to  10  KHz  was  examined  comparing  the  different  interpolation  functions.  A  fre¬ 
quency  sweep  step  size  of  50  Hz  was  chosen,  and  the  relative  error  is  shown  in  Figure  8.2.2. 
While  the  4^'  Lagrangian  and  the  2'’*  Hermitian  interpolation  functions  show  relatively  large  er¬ 
rors,  the  dP'  Hermitian  interpolation  function  has  very  low  error.  The  comparison  of  the  %  error 
of  the  magnitude  of  the  displacement  for  the  “master”  DOF  u„  is  shown  in  Figure  8.2.3.  The  4”' 
Lagrangian  and  the  2^  Hermitian  interpolation  methods  show  very  large  errors  (approximately 
10%)  compared  to  the  full  DOF  solution  in  Figure  8.2.4.  Where  as  the  4'”  Hermitian  interpolation 
shows  a  surprisingly  low  error  of  approximately  0.01%.  All  the  interpolation  methods  seem  to 
capture  the  eigenvalues  that  are  found  within  the  frequency  window,  Table  8.2.1. 

Table  8.2.1  -  Calculated  Constrained  Eigenvalues 

0.2652 

2.0918 

2.1092 

4.6078 

11.6655 

15.2256 

16.4790 

18.9294 

22.0093 

23.1288 
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From  both  Figures  8.2.2  and  8.2.3,  it  can  be  seen  that  errors  are  higher  in  the  low  fre¬ 
quency  range  and  decrease  across  the  frequency  window.  This  fact  arises  from  the  presence  of  an 
eigenvalue  at  the  very  low  end  of  the  window.  The  presence  of  this  eigenvalue  has  an  effect  on 
the  smoothness  of  the  interpolated  function  in  the  low  end  of  the  frequency  window.  This  shows 
that  the  interpolated  function  is  sensitive  to  the  eigenvalues  of  the  system  particularly  if  they  are 
at  or  near  the  end  interpolation  (inversion)  points.  Several  options  are  available  to  improve  the 
over  all  accuracy;  shifting  of  the  frequency  interpolation  window  away  from  the  end  point  eigen¬ 
values,  the  reduction  of  the  frequency  window  range  or  variable  placement  the  interior  interpola¬ 
tion  (inversion)  points. 

It  is  possible  within  the  FWR  codes  to  vary  the  placement  of  the  interior  interpolation 
(inversion)  points  for  both  the  Lagrangian  and  Hermitian  interpolation  methods.  The  above  fre¬ 
quency  window  was  examined  by  placing  the  interior  points  at  1.5  KHz  and  7.5  KHz  (Figures 
8.2.5  to  8.2.7)  and  at  1.5  KHz  and  4  KHz  (Figures  8.2.8  to  8.2.10).  While  the  movement  of  the 
interior  points  does  not  generally  reduce  the  overall  error,  it  does  however  tend  to  smooth  the 
error  more  uniformly  across  the  entire  window.  This  can  most  readily  be  seen  in  the  %  error  of 
the  magnitude  of  the  displacement  for  the  “master”  DOF  u,,  in  Figures  8.2.6  and  8.2.9  as  com¬ 
pared  to  Figure  8.2.3  (interior  points  at  3.3  KHz  and  6.7  KHz).  It  can  be  seen  in  all  three  figures, 
that  the  %  error  seems  to  be  apparently  the  same  at  1%  and  0.001%  for  the  Lagrangian  and  for 
the  dP'  Hermitian  interpolation  functions  respectively.  However  it  is  clear  that  the  variable 
placement  of  the  interior  point  does  effect  the  error  in  different  smaller  local  frequency  regions 
within  the  window.  There  is  apparently  no  significant  effect  on  the  actual  solution  and  its  accu¬ 
racy,  i.e.  magnitude  of  the  displacement  for  a  particular  node  point. 

A  further  reduced  frequency  window  range,  0  to  5  KHz,  was  examined  for  a  frequency 
sweep  step  size  of  25  Hz.  This  reduced  window  both  improves  the  over  all  accuracy  as  well  as 
reduces  the  interpolation  error  associated  with  the  eigenvalue  at  the  lower  frequency  interpola¬ 
tion  point.  Figure  8.2.11.  The  slightly  higher  error  at  the  low  frequencies  is  still  present  but  is 
clearly  reduced.  The  %  error  of  the  magnitude  of  the  displacement  of  DOF  shown  in  Figure 
8.2.12  also  has  reduced  error  for  both  the  d^'  Lagrangian  and  the  2'”  Hermitian  interpolation 
functions  to  approximately  1%  and  d*”  Hermitian  interpolation  function  has  a  much  lower  error  at 
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0.00001%.  As  seen  in  Figure  8.2.13,  the  magnitude  of  the  displacement  of  DOF  u,,,  now 
compares  extremely  well  to  the  full  DOF  solution  for  both  the  4^'  Lagrangian  and  the  4^' 
Hermitian  interpolation  functions. 
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3-Layer  Steel-Polystyrene  2D  Figures 


Figure  8.2.1  -  3-Layer  Steel-Polystyrene  Laminate  2D  Model 


Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 


Figure  8.2.2  -  Relative  Error  comparing  different  interpolation  functions  for  the  frequency  range  of  0  to 
10  KHz  and  a  frequency  step  size  of  50  Hz. 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
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Figure  8.2.3  -  %  Error  for  the  “master”  DOF  luli^-comparing  interpolations  for  the  frequency  range  of  0  to  10 
KHz  and  a  step  size  of  50  Hz.  NOTE:  Interior  interpolation  (inversion)  points  are  at  3.3  KHz  and  6.7  KHz. 


•  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  Displacement  Magnitude 


Figure  8.2.4  -  The  “master”  DOF  lulj^  (magnitude  of  the  displacement)  comparing  two  interpolation 
functions  to  the  full  DOF  solution  for  the  frequency  range  of  0  to  10  KHz. 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  Displacement  Magnitude  %  Error 
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Figure  8.2.6  -  %  Error  for  the  “master”  DOF  lulj^-comparing  interpolations  for  the  range  of  0  to  10  KHz  and  a  step 
size  of  50  Hz.  NOTE:  Interior  interpolation  (inversion)  points  are  placed  at  1.5  KHz  and  7.5  KHz. 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  Displacement  Magnitude 
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Figure  8.2.7  -  The  “master”  DOF  lul,,  (magnitude  of  the  displacement)  comparing  the  two  interpolations  to 
the  full  DOF  solution.  NOTE:  Interior  interpolation  (inversion)  points  are  placed  at  1.5  KHz  and  7.5  KHz. 


ayer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 


Figure  8.2.8  -  Relative  Error  comparing  interpolation  functions  for  the  range  of  0  to  10  KHz  and  a  frequency 
step  size  of  50  Hz.  NOTE:  Interior  interpolation  (inversion)  points  are  placed  at  1.5  KHz  and  4  KHz. 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  u  Displacement  Magnitude  %  Error 
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Figure  8.2.9  -  %  Error  for  the  “master”  DOF  luli^-comparing  interpolation  functions  for  the  range  of  0  to  10 
laiz  and  a  step  size  of  50  Hz.  NOTE:  Interior  interpolation  (inversion)  points  are  at  1.5  KHz  and  4  KHz. 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  u  Displacement  Magnitude 


Figure  8.2.10  -  The  “master”  DOF  lulj^  (magnitude  of  the  displacement)  comparing  the  two  interpolations 
to  the  full  DOF  solution.  NOTE:  Interior  interpolation  (inversion)  points  are  placed  at  1 .5  KHz  and  4  KHz. 


^ayer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  Displacement  Magnitude  %  Error 


Figure  8.2.12  -  %  Error  for  the  “master”  DOF  lul,^  comparing  interpolations  for  the  frequency  range  of  0 
to  5  KHz  and  a  frequency  step  size  of  25  Hz. 


3  Layer  Laminate  :  230  DOF  to  70  DOF  (70%  Reduction) 
Master  DOF  Uj^^  Displacement  Magnitude 


Figure  8.2.13  -  The  “master”  DOF  lul,^  (magnitude  of  the  displacement)  comparing  the  4pt  Lagrangian  and 
the  4pt  Hermitian  interpolation  function  to  the  full  DOF  solution  for  the  frequency  range  of  0  to  5  KHz. 


8.3  Spherical  Half-Shell  Model 


The  next  problem  considered  is  a  simple  spherical  shell  as  shown  in  Figure  8.3.1  with  the 
properties  listed  in  Table  8.3.1.  A  similar  problem  is  discussed  in  the  ABAQUS™  User's  Manual 
[10].  The  model  is  constructed  using  50  quadratic  axisymmetric  SAX2  isoparametric  shell  ele¬ 
ments  for  the  undamped,  free  vibration  problem.  The  total  number  of  node  points  is  101  with  3 
degrees  of  freedom  per  node  point;  hence  the  total  DOFs  for  this  model  is  303.  The  number  of 
“master”  node  points  was  chosen  as  29,  so  the  degree  of  reduction  of  the  model  is  71%.  The  se¬ 
lected  “master”  node  points  are  shown  in  the  Figure  8.3.1  as  “•”  symbols  (Node  #s  1,  51  and 

101)  which  represent  the  loaded  and  constrained  nodes.  The  remaining  “master”  nodes  were  ar¬ 
bitrarily  chosen  from  the  model  (Node  #s  5, 9,  13, 17, 21, 25, 27, 29,  33  and  37). 

Table  8.3.1  -  Steel  Half-Spherical  Shell  Properties 


Young’s  Modulus 

180.0  GPa 

Poisson  Ratio 

0.333 

Density 

7670  Kg/m' 

Initially  a  large  frequency  interpolation  window  of  0  to  500  Hz  was  examined  with  a 
fairly  coarse  frequency  sweep  step  size  of  2  Hz.  Only  the  d*”*  Lagrangian  and  the  4'”  Hermitian 
interpolation  methods  were  compared,  shown  in  Figure  8.3.2.  As  can  be  seen  there  are  a  large 
number  of  eigenvalues  in  this  window.  Surprisingly  both  interpolation  functions  capture  the 
functional  behavior  at  all  of  the  eigenvalues  despite  the  large  window  and  coarse  step  size.  In 
Figure  8.3.3  the  Lagrangian  interpolation  function  shows  a  larger  %  error  of  the  magnitude  of  the 
displacement  of  the  “master”  node  Usi^  of  approximately  10%,  where  as  the  Hermitian  shows  an 
approximate  error  of  only  0.01%.  From  this  figure  and  the  actual  magnitude  of  the  displacement 
of  the  “master”  node  Uj,,,  of  Figure  8.3.4,  it  can  clearly  be  seen  that  the  functional  behavior  at 
every  eigenvalue  is  represented.  The  Figures  8.3.5  and  8.3.6  show  the  %  error  of  the  magnitude 
of  the  displacement  of  another  “master”  node  Ujoiy  and  its  magnitude.  Identical  errors  as  well  as 
the  direct  parallel  response  at  the  eigenvalues  are  seen  as  above.  Figures  8.3.7  to  8.3.12  show  the 
%  error  of  the  magnitude  of  the  displacement  and  the  magnitude  of  the  displacements  for  a  “re¬ 
constructed”  “slave”  node  U75  for  each  of  its  three  degrees  of  freedom.  All  of  these  plots  demon- 
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strate  the  consistency  of  the  error  and  the  accuracy  of  the  response  of  the  reconstructed  solution 
as  compared  to  the  full  solution  for  this  model. 

The  interpolation  frequency  window  size  was  reduced  to  0  to  200  Hz  with  the  sweep  step 
size  to  1  Hz.  Now  within  this  particular  window  only  two  eigenvalues  are  present  near  the  end 
points  of  the  window.  The  relative  error,  Figure  8.3.13  and  the  %  error  of  the  “master”  node  Uj,^, 
Figure  8.3.14,  show  a  dramatic  decrease  and  very  accurate  behavior  at  the  eigenvalues.  The  4'" 
Hermitian  seems  to  have  extreme  accuracy  in  this  frequency  range.  The  actual  magnitude  of  the 
displacements  for  this  DOF  shows  no  difference  in  the  reduced  interpolated  solutions  and  the  full 
DOF  solution.  Figure  8.3.15.  These  figures  were  for  a  level  of  71%  DOF  reduction  (303  DOFs  to 
87  DOFs).  This  same  window  range  was  then  examined  with  a  higher  level  of  DOF  reduction  of 
303  DOFs  to  29  DOFs  or  an  87%  DOF  reduction.  Surprisingly  the  relative  errors.  Figure  8.3.16 
and  the  %  errors.  Figure  8.3.17,  change  very  little,  increasing  by  approximately  one  order  of 
magnitude.  Careful  examination  of  these  two  figures  and  the  magnitude  of  the  displacement  in 
Figure  8.3.18  shows  only  that  the  4'*'  Lagrangian  does  not  match  exactly  the  eigenvalue  behavior 
at  the  low  range  of  the  frequency  window.  However,  the  levels  of  error  even  for  this  significant 
DOF  reduction  are  quite  surprising  at  1%  for  the  Lagrangian  and  0.0001%  for  the  Hermitian 
methods. 

As  is  seen  in  Figure  8.3.2  or  Figure  8.3.19,  there  are  a  large  number  of  eigenvalues  within 
the  frequency  range  of  200  Hz  to  300  Hz.  The  frequency  window  range  of  200  to  300  Hz  was 
examined  at  both  the  71%  DOF  reduction  and  the  87%  DOF  reduction  level.  This  window  range 
was  chosen  in  order  to  characterize  the  behavior  of  the  FWR  method  with  a  large  number  of 
eigenvalues  with  in  a  small  frequency  window.  Figure  8.3.20  shows  the  relative  error  at  71%  re¬ 
duction  for  the  window  range  and  the  sweep  step  size  of  0.5  Hz.  Both  interpolation  methods  re¬ 
sult  in  very  low  %  errors  («0.1%)  and  highly  accurate  responses  as  compare  to  the  full  DOF  so¬ 
lution,  Figures  8.3.21  and  8.3.22.  The  characterization  at  87%  DOF  reduction  is  shown  in 
Figures  8.3.23  to  8.3.25.  These  figures  show  identical  trends  as  described  above.  It  should  be 
noted  that  as  in  the  above  cases  discussed,  all  eigenvalue  responses  are  clearly  captured  in  both 
DOF  reductions. 
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An  important  fact  can  be  observed  with  respect  to  the  4”'  Hermitian  interpolation  curves 
in  these  figures.  This  is  that  the  end  interpolation  (inversion)  points  of  the  frequency  window  do 
not  approach  a  nearly  zero  value,  as  does  the  Lagrangian  interpolation.  It  is  important  to  note  that 
the  Hermitian  interpolation  method  includes  derivative  information  at  all  of  the  interpolation 
points.  As  a  result  if  an  eigenvalue  is  near  the  interpolation  point  particularly  at  the  frequency 
window  end  points,  the  function  and  its  derivative  are  rapidly  varying  near  the  interpolation 
points  in  region  of  validity  for  the  interpolation  function.  For  this  window  range  of  200  Hz  to 
300  Hz,  it  can  be  seen  from  Table  8.3.2  (or  Figure  8.3.19)  that  there  are  eigenvalues  that  lie  just 
outside  of  the  interpolation  window  range  (188.14  Hz  and  308.71  Hz).  Hence  the  influence  from 
these  eigenvalues,  even  though  outside  of  the  interpolation  region,  does  in  fact  contribute  to  the 
smoothness  of  the  interpolation  function  and,  hence,  to  the  accuracy  of  the  response.  This  shows 
that  the  continuity  of  the  interpolation  function  will  be  maintained  across  the  interpolation  win¬ 
dow  end  points.  The  Lagrangian  value  should  always  be  nearly  zero  at  an  interpolation  point. 


Table  8.3.2  -  Calculated  Eigenvalues 


iSi;Wnconstrained 

187.36  Hz 

188.14  Hz 

222.69  Hz 

223.48  Hz 

236.95  Hz 

237.73  Hz 

244.41  Hz* 

245.23  Hz 

249.30  Hz 

250.23  Hz 

253.30  Hz* 

254.44  Hz 

257.26  Hz 

258.75  Hz 

267.02  Hz 

263.68  Hz 

273.53  Hz 

269.63  Hz 

281.49  Hz 

276.87  Hz 

291.12  Hz 

285.68  Hz 

302.59  Hz 

296.24  Hz 

316.06  Hz 

308.71  Hz 

*  Used  for  Rayleigh  Damping  Calculation 


Another  very  important  note  that  applies  to  all  previously  described  examples  is  that  the 
eigenvalues  of  the  system  were  NOT  known  beforehand.  Having  shown  how  accurately  the  re- 
sponse  of  the  reduced  DOF  problem  can  be  obtained,  it  is  clear  that  the  eigenvalues  of  the  full 
DOF  system  may  be  closely  approximated  from  the  reduced  DOF  model  analysis.  It  is  important 
to  note  that  these  eigenvalues  represent  the  resonances  of  the  “reduced”  DOF  system,  which 
generally  have  been  found  to  be  not  significantly  different  from  those  of  the  full  DOF  system.  In 


90 


order  to  obtain  more  accurate  eigenvalues  of  the  reduced  system,  small  frequency  sweep  step 
sizes  would  be  necessary  which  could  a  performed  by  a  partial  frequency  sweep  within  an 
interpolated  window.  As  described  above,  this  is  where  a  window  range  is  interpolated  and  then 
only  a  small  portion  of  that  interpolation  window  is  swept  at  a  very  small  step  size.  (See  for 
example  Figures  2.2. 1.6  and  2.2. 1.7)  These  features  could  be  extremely  useful  in  developing 
models,  understanding  system  responses,  obtaining  eigenvalues  and,  most  importantly,  in 
guiding  the  analysis  of  large  complex  structural  models. 

One  additional  aspect  of  FWR  method  is  that  multiple  frequency  windows  may  be  ana¬ 
lyzed  and  these  solutions  can  be  then  combined  to  construct  a  larger  frequency  response  window. 
To  illustrate  this  feature  of  FWR  method,  a  frequency  interpolation  window  of  0  to  300  Hz  was 
analyzed  with  a  frequency  sweep  step  size  of  1  Hz  at  71%  DOF  reduction.  These  results  are  pre¬ 
sented  in  Figures  8.3.26  to  8.3.28.  The  above  analyses  for  the  frequency  windows  of  0  to  200  Hz 
and  200  to  300  Hz,  can  be  combined  and  compared  to  the  full  analysis  for  the  window  range  of  0 
to  300  Hz.  The  combined  analyses  for  the  window  range  of  0  to  200  Hz,  Figures  8.3.13  to  8.3.15, 
with  that  for  200  to  300  Hz  window  range.  Figures  8.3.20  to  8.3.22  are  shown  in  Figures  8.3.29 
to  8.3.3 1 .  It  is  difficult  to  compare  the  overall  accuracy  (error)  of  the  combined  data  as  compared 
to  the  full  frequency  range  because  the  two  analyses,  which  are  combined,  have  different  pa¬ 
rameters  (window  and  step  size).  It  is  clear  that  the  responses  of  the  magnitude  of  the  displace¬ 
ments  shown  in  Figure  8.3.28  (full  window)  and  Figure  8.3.31  (combined  windows)  are  identi¬ 
cal.  The  combinations  of  the  0  to  200  Hz  window  for  87%  DOF  reduction  (Figures  8.3.16  to 
8.3.18),  with  that  for  200  to  300  Hz  window  range  at  71%  reduction  (Figures  8.3.20  to  8.3.22) 
are  shown  in  Figures  8.3.32  to  8.3.34.  These  figures  also  show  nearly  identical  response  as  com¬ 
pare  to  the  full  0  to  300  Hz  window,  despite  representing  different  degrees  of  reduction.  Finally  a 
comparison  of  all  the  responses  of  the  magnitude  of  the  displacements  (Figures  8.3.28,  8.3.31 
and  8.3.34)  to  the  analysis  of  the  original  0  to  500  Hz  window.  Figure  8.3.35,  also  demonstrates 
the  consistency  of  this  method.  Note  that  the  scale  for  Figure  8.3.35  has  been  redrawn. 

Rayleigh  damping  was  next  introduced  into  the  Spherical  Half  Shell  Model.  The  mass 
and  stiffness  matrices  were  obtained  from  an  ABAQUS™  model  and  reordered  according  to  the 
selection  of  the  “master”  nodes.  These  matrices  were  then  used  to  calculate  the  Rayleigh  damp- 
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ing  as  defined  as 


D  s  aM  +  pS 
Where  the  terms 


(8.3.1) 


P  _  ~  ^1^1 )  (8.3.2) 

And 

a  =  2o}^^^  -  pcOi  =  to,  (2^,  -  )3<0j)  (8.3.3) 


Here  the  terms,  ^^nd  4  are  the  weightings  or  percentages  of  critical  damping  of  the  select  fre¬ 
quency  modes,  tOi  and  <W2.  From  the  calculated  eigenvalues  of  this  model  listed  in  Table  8.3.2, 
the  frequency  modes  of  o)^  =  244.41  Hz  and  0)^  =  253.30  Hz  were  chosen  as  being  approxi¬ 
mately  in  the  middle  of  the  frequency  window  range  of  200  Hz  to  300  Hz.  The  calculated  values. 
Table  8.3.3  for  a  and  p,  according  to  equations  8.3.2  and  8.3.3,  at  these  selected  eigenvalues 
are  summarized  for  different  levels  of  critical  damping  (5%,  1%  and  0.5%). 


Table  8.3.3  -  Calculated  Values  for  a  and  p  at  selected  eigenvalues 
m,  =  244.41  Hz  and  0)2  =  253.30  Hz 


4 

4 

a 

P 

0.05 

0.05 

78.1552 

3.1977  X  lO'® 

0.01 

0.01 

15.6310 

6.3955  X  10  " 

0.005 

0.005 

7.81552 

3.1977  X  10  " 

Initially  5%  of  critical  damping  was  examined  for  this  model  for  the  frequency  window 
range  of  200  Hz  to  300  Hz  and  a  frequency  sweep  step  of  0.25  Hz.  It  was  found,  as  expected, 
that  both  real  and  imaginary  solution  responses  are  obtained.  Plots  of  the  relative  error  for  the 
real  and  imaginary  solutions  as  compared  to  the  full  solution  show  very  small  relative  error. 
Figures  8.3.36  and  8.3.37.  It  is  important  to  note  that  the  response  near  the  eigenvalues  does  not 
appear  as  in  the  case  with  no  damping  (Figure  8.3.20).  Both  the  real  and  imaginary  %  error  for 
the  magnitude  of  the  displacement  of  DOF  U5i,t  are  shown  in  Figures  8.3.38  and  8.3.39.  These 
figures  show  an  error  of  less  than  0.01%.  Comparison  of  this  %  error  in  the  case  in  which 
damping  is  present  to  the  case  of  no  damping  (Figures  8.3.20  to  8.3.22)  show  near  the  same  % 


92 


error.  Comparison  of  the  real  and  imaginary  responses  of  the  magnitude  of  the  displacement  of 
the  “master”  DOF  seen  in  Figures  8.3.40  and  8.3.41  shows  that  the  reduced  methods  as  in¬ 
distinguishable  from  the  full  solution  response.  Again  most  notable  in  all  of  these  figures  is  the 
apparent  lack  of  influence  from  the  eigenvalues.  This  observation  produced  some  concern  as  it 
was  assumed  that  a  small  fraction  of  critical  damping  (5%)  should  only  tend  to  smooth  and  lower 
the  peaks  associated  with  eigenvalues. 


Table  8.3.4  -  Calculated  Constrained  Eigenvalues 


No  Damping 

11.2585438  Hz 

9.80502865  Hz 

188.137713  Hz 

188.098833  Hz 

223.481103  Hz 

223.34894  Hz 

237.734345  Hz 

>237.547289flz 

245.229167  Hz 

245.008992'Hz 

250.226527  Hz 

249.982589  Hz 

254.441903,  Hz 

254.176854  Hz  < .  • 

258.748788  Hz  . 

258.461147  Hz  ' 

i*;:,.  263.68^2892  Hz,:,:  rfj 

|Pv285.(^79'29l: 
lf#?'296.240074  Hzffi^ 

.  ,>263368079  Hz 

^§76i480859ilz 

,  -■'i^225972:H|:f  *' 

.  29Xt09751  Hz 

308.719867  Hz 

308.08933  Hz 

323.24202  Hz 

322.481799  Hz 

339.893522  Hz 

338.966522  Hz 

358.728968  Hz 

357.588881  Hz 

379.775347  Hz 

378.364823  Hz 

403.035042  Hz 

401.283648  Hz 

428.459824  Hz 

426.282331  Hz 

445.0769  Hz 

442.589095  Hz 

456.411232  Hz 

453.696729  Hz 

486.150759  Hz 

482.780906  Hz 

In  order  to  try  to  clarify  the  influence  of  the  damping  as  related  to  the  eigenvalue  re¬ 
sponse,  a  smaller  percentage  of  the  critical  damping  of  1%  was  chosen.  As  is  seen  in  Figure 
8.3.42  for  the  real  relative  error  and  in  Figure  8.3.43  for  the  imaginary  relative  error,  the  re¬ 
sponse  near  the  eigenvalues  is  clearly  shown.  The  comparison  of  Figure  8.3.42  to  Figure  8.3.20, 
the  relative  error  of  the  case  with  no  damping,  as  expected,  clearly  indicates  that  the  peaks  near 
the  eigenvalues  are  reduced  in  height  and  more  rounded.  Hence,  at  1%  of  the  critical  damping 
the  obvious  influence  of  the  eigenvalues  is  present.  The  %  error  for  both  the  real  and  imaginary 
magnitude  of  the  displacement  of  the  master  DOF  is  shown  in  Figures  8.3.44  and  8.3.45.  It  is 
interesting  to  note  that  the  %  error  of  «0.01%  is  nearly  identical  to  that  of  the  undamped  case 
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(Figure  8.3.21).  Also,  as  in  the  undamped  case  (Figure  8.3.22),  there  is  no  discemable  difference 
between  the  reduced  response  and  the  full  solution  response  in  Figures  8.3.46  and  8.3.47.  The 
further  reduction  in  the  percentage  of  the  critical  damping  of  0.5%  was  examined.  Figures  8.3.48 
to  8.3.53  show  identical  results  to  the  previous  damping  cases.  The  peak  heights  and  peak 
rounding  are  less  than  the  undamped  case  but  slightly  more  than  the  1%  damping  case. 

The  comparison  of  the  most  stringent  error  criteria,  the  relative  error,  for  the  case  in 
which  there  is  no  damping  to  the  cases  with  different  percentages  of  critical  damping  demon¬ 
strates  the  effect  that  damping  has  on  the  overall  response  of  this  model  particularly  the  response 
near  the  eigenvalues.  The  relative  error  for  the  undamped  case  is  compared  to  that  of  the  three 
damped  cases  (real  component  only)  for  the  Lagrangian  (Figure  8.3.54)  and  the  Hermitian  inter¬ 
polation  methods  (Figure  8.3.55).  These  figures  show  clearly  the  gradual  decrease  and  smooth¬ 
ing  if  the  peaks  at  (near)  the  eigenvalues.  Also  the  comparison  of  the  %  error  of  the  magnitude  of 
the  displacements  of  the  “master”  DOF  U5,,,  for  the  undamped  and  damped  cases.  Figures  8.3.21, 
8.3.28,  8.3.44  and  8.3.50,  shows  that  the  overall  error  is  the  same  at  0.01%.  It  appears  that  the 
accuracy  of  the  FWR  method  is  related  to  the  %  DOF  reduction,  the  interpolation  window  size, 
the  sweep  step  size  and  interpolation  method  used  and  not  whether  damping  is  present.  This 
could  allow  analyses  to  be  performed  with  out  damping  in  order  to  determine  the  important  pa¬ 
rameters  for  the  analysis  and  then  used  to  examine  cases  with  damping. 

It  is  recognized  that  perhaps  more  appropriate  choices  of  the  frequency  modes  used  to 
calculate  the  Rayleigh  damping  for  tiiis  model  could  be  made.  Initially  the  frequency  modes  cho¬ 
sen  were  approximately  in  the  middle  of  the  interpolation  window.  These  were  selected  not  nec¬ 
essarily  to  represent  a  real  problem  with  damping  but  rather  these  were  preliminary  choices  used 
to  determine  the  applicability  of  the  FWR  method  to  systems  with  Rayleigh  damping. 
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Spherical  Half-Shell  Model  Figures 


Figure  8.3.1  -  Half-Spherical  Shell  Model 


Half  Sphere :  303  DOF  to  87  DOF  (71%  Reduction) 


Figure  8.3.2  -  Relative  Error  comparing  interpolations  for  the  frequency  range  of  0  to  500  Hz  and  a 
frequency  step  size  of  2  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Master  DOF  u Displacement  Magnitude  %  Error 
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Figure  8.3.5  -  %  Error  for  the  “master”  DOF  lul,oiy  comparing  interpolations  for  the  range  of  0  to  500  Hz 
and  a  frequency  step  size  of  2  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Master  DOF  Ui  Displacement  Magnitude 
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Figure  8.3.6  -  The  “master”  DOF  lul,oiy  (magnitude  of  the  displacement)  comparing  interpolation 
functions  for  the  frequency  range  of  0  to  500  Hz  and  a  frequency  step  size  of  2  Hz. 


Half  Sphere :  303  DOF  to  87  DOF  (71%  Reduction) 
Slave  DOF  u,^^^  Displacement  Magnitude  %  Error 
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Figure  8.3.7  -  %  Error  for  the  reconstructed  “slave”  DOF  IUI75X  comparing  interpolations  for  the  range  of  0 
to  500  Hz  and  a  step  size  of  2  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Slave  DOF  Displacement  Magnitude 
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Figure  8.3.8  -  The  reconstructed  “slave”  DOF  lu^j^  for  the  frequency  range  of  0  to  500  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Slave  DOF  07^  Displacement  Magnitude 


Figure  8.3.10  -  The  reconstructed  “slave”  DOF  lu^  for  the  frequency  range  of  0  to  500  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Slave  DOF  U75j.  Displacement  Magnitude 
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Figure  8.3.12  -  The  reconstructed  “slave”  DOF  lu^j,  for  the  frequency  range  of  0  to  500  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
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Figure  8.3. 13  -  Relative  Error  comparing  interpolations  for  the  range  of  0  to  200  Hz  and  a  step  size  of  1  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Master  DOF  Displacement  Magnitude 
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Figure  8.3.15  -  The  “master”  DOF  lulj,,  (magnitude  of  the  displacement)  comparing  interpolation 
functions  for  the  range  of  0  to  200  Hz. 


Half  Sphere :  303  DOF  to  39  DOF  (87%  Reduction) 


Figure  8.3.16  -  Relative  Error  comparing  interpolation  functions  for  the  frequency  range  of  0  to  200  Hz 
and  a  frequency  step  size  of  1  Hz.  NOTE:  87%  DOF  Reduction. 


Half  Sphere  :  303  DOF  Constrained  System  Eigenevalues 
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Figure  8.3.19  -  The  constrained  eigenvalues  for  the  spherical  half  shell  model. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
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Figure  8.3.20  -  Relative  Error  comparing  interpolation  functions  for  the  frequency  range  of  200  to  300  Hz 
and  a  frequency  step  size  of  0.5  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Master  DOF  u^^x  Displacement  Magnitude  %  Error 
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Figure  8.3.21  -  %  Error  for  the  “master”  DOF  luls,^  comparing  interpolation  functions  for  the  frequency  range 
of  200  to  300  Hz  and  a  frequency  step  size  of  0.5  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Master  DOF  Displacement  Magnitude 


Figure  8.3.22  -  The  “master”  DOF  lulji,  comparing  interpolation  functions  for  the  frequency  range  of  200 
to  300  Hz  and  a  frequency  step  size  of  0.5  Hz. 


Half  Sphere  :  303  DOF  to  39  DOF  (87%  Reduction) 
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Figure  8.3.23  -  Relative  Error  comparing  interpolation  functions  for  the  frequency  range  of  200  to  300  Hz 
and  a  frequency  step  size  of  0.5  Hz.  NOTE;  87%  DOF  Reduction. 


Half  Sphere  :  303  DOF  to  39  DOF  (87%  Reduction) 
Master  DOF  Ug  Displacement  Magnitude 
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Figure  8.3.25  The  “master”  DOF  (magnitude  of  the  displacement).  NOTE:  87%  DOF  Reduction. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction)  0:300  Hz 
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Figure  8.3.26  -  Relative  Error  comparing  different  interpolation  functions  for  the  frequency  range  of  0  to 
300  Hz  and  a  frequency  step  size  of  1  Hz. 


Half  Sphere :  303  DOF  to  87  DOF  (71%  Reduction)  0:300  Hz 
Master  DOF  Displacement  Magnitude 


Figure  8.3.28  The  “master”  DOF  lulji^  for  the  frequency  range  of  0  to  300  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 


Figure  8.3.29  -  Relative  Error  of  combined  data  0  to  200  Hz  and  200  to  300  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Master  DOF  u^^^  Displacement  Magnitude  %  Error 
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Figure  8.3.30  -  %  Error  of  combined  data  0  to  200  Hz  and  200  to  300  Hz  at  71%  DOF  reduction. 


Half  Sphere :  303  DOF  to  39  DOF  (87  %  Reduction)  0  to  200  Hz 
Half  Sphere  :  303  DOF  to  87  DOF  (71  %  Reduction)  200  to  300  Hz 


Figure  8.3.32  -  Relative  Error  of  combined  data  0  to  200  Hz  at  87%  DOF  reduction 
71%  DOF  reduction. 


Half  Sphere :  303  DOF  to  39  DOF  (87%  Reduction)  0  to  200  Hz 
Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction)  200  to  300  Hz 
Master  DOF  u^  Displacement  Magnitude 


to  300  Hz  at  71%  DOF  reduction. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction)  0:500  Hz 
Master  DOF  Displacement  Magnitude  %  Error 


Figure  8.3.35  -  The  “master”  DOF  lulj,,  for  the  data  0  to  500  Hz  at  71%  DOF  reduction.  NOTE:  Expanded 
frequency  scale  from  Figure  8.3.4. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
5%  Rayleigh  Damping  -  Real  Component 
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Figure  8.3.36  -  Relative  Error,  real  component,  comparing  interpolations  for  the  frequency  range  of  200 
to  300  Hz  and  a  frequency  step  size  of  0.5  Hz.  NOTE:  5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
5%  Rayleigh  Damping  -  Imaginary  Component 


Figure  8.3.37  -  Relative  Error,  imaginary  component,  comparing  interpolations  for  the  frequency  range  of 
200  to  300  Hz  and  a  frequency  step  size  of  0.5  Hz.  NOTE:  5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
5%  Rayleigh  Damping  -  Imaginary  Component 
Master  DOF  u^ Displacement  Magnitude  %  Error 
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Figure  8.3.39  -  %  Error,  imaginary  component,  for  the  “master”  DOF  lulj,^  comparing  interpolations  for  the 
frequency  range  of  200  to  300  Hz.  NOTE:  5%  of  critical  damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
5%  Rayleigh  Damping  -  Real  Component 
Master  DOF  Displacement  Magnitude 
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size  of  0.5  Hz.  NOTE:  5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
5%  Rayleigh  Damping  -  Imaginary  Component 
Master  DOF  Displacement  Magnitude 
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Figure  8.3.41  -  The  imaginary  component  of  the  “master”  DOF  for  the  frequency  range  of  200  to  300  Hz 
and  a  frequency  step  size  of  0.5  Hz.  NOTE:  5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
1%  Rayleigh  Damping  •  Real  Component 


Figure  8.3.42  -  Relative  Error,  real  component,  comparing  interpolations  for  the  frequency  range  of  200 
to  300  Hz  and  a  frequency  step  size  of  0.25  Hz.  NOTE;  1%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
1%  Rayleigh  Damping  ■  Imaginary  Component 
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Figure  8.3.43  -  Relative  Error,  imaginary  component,  comparing  interpolations  for  the  frequency  range 
of  200  to  300  Hz  and  a  frequency  step  size  of  0.25  Hz.  NOTE:  1%  of  critical  damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
1%  Rayleigh  Damping  -  Imaginary  Component 
Master  DOF  Ug Displacement  Magnitude  %  Error 
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Figure  8.3.45  -  %  Error,  imaginary  component,  for  the  “master”  DOF  lulj,,,  comparing  interpolations  for 
the  frequency  range  of  200  to  300  Hz.  NOTE:  1%  Damping. 


Half  Sphere :  303  DOF  to  87  DOF  (71%  Reduction) 
1%  Rayleigh  Damping  -  Real  Component 
Master  DOF  Displacement  Magnitude 


Figure  8.3.46  -  The  real  component  of  the  “master”  DOF  lulj,,  for  the  frequency  range  of  200  to  300  Hz 
and  a  frequency  step  size  of  0.25  Hz.  NOTE:  1%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
1%  Rayleigh  Damping  -  Imaginary  Component 
Master  DOF  Displacement  Magnitude 
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Figure  8.3.47  -  The  imaginary  component  of  the  “master”  DOF  lulji^  for  the  frequency  range  of  200  to 
300  Hz  and  a  frequency  step  size  of  0.25  Hz.  NOTE:  1%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
0.5%  Rayleigh  Damping  -  Real  Component 
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Figure  8.3.48  -  Relative  Error,  real  component,  comparing  interpolation  functions  for  the  frequency  range 
of  200  to  300  Hz  and  a  frequency  step  size  of  0.25  Hz.  NOTE;  0.5%  of  critical  damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
0.5%  Rayleigh  Damping  -  Imaginary  Component 


Figure  8.3.49  -  Relative  Error,  imaginary  component,  comparing  interpolation  functions  for  the  frequency 
range  of  200  to  300  Hz  and  a  frequency  step  size  of  0.25  Hz.  NOTE:  0.5%  Damping. 


Half  Sphere :  303  DOF  to  87  DOF  (71%  Reduction) 
0.5%  Rayleigh  Damping  -  Real  Component 
Master  DOF  Displacement  Magnitude  %  Error 
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Figure  8.3.50  -  %  Error,  real  component,  for  the  “master”  DOF  lulj,,^  comparing  interpolations  for  the 
frequency  range  of  200  to  300  Hz.  NOTE:  0.5%  Damping  with  less  than  0.01%  Error. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
0.5%  Rayleigh  Damping  -  Imaginary  Component 
Master  DOF  Displacement  Magnitude  %  Error 
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Figure  8.3.51  -  %  Error,  imaginary  component,  for  the  “master”  DOF  comparing  interpolations  for 
the  frequency  range  of  200  to  300  Hz.  NOTE:  0.5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
0.5%  Rayleigh  Damping  -  Real  Component 
Master  DOF  Displacement  Magnitude 


Figure  8.3.52  -  The  real  component  of  the  “master”  DOF  lulj,^  for  the  frequency  range  of  200  to  300  Hz 
and  a  frequency  step  size  of  0.25  Hz.  NOTE:  0.5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
0.5%  Rayleigh  Damping  -  Imaginary  Component 
Master  DOF  Displacement  Magnitude 
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Figure  8.3.53  -  The  imaginary  component  of  the  “master”  DOF  lulj,^.  NOTE:  0.5%  Damping. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Comparison  of  Lagrangian  Method  for  Different  Damping 


Figure  8.3.54  -  Relative  Error,  real  component,  comparing  level  of  critical  damping  for  Lagrangian 
interpolation  method  for  the  frequency  range  of  200  to  300  Hz  and  a  frequency  step  size  of  0.5  Hz. 


Half  Sphere  :  303  DOF  to  87  DOF  (71%  Reduction) 
Comparison  of  Hermitian  Method  for  Different  Damping 
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Figure  8.3.55  -  Relative  Error,  real  component,  comparing  level  of  critical  damping  for  Hermitian 
interpolation  method  for  the  frequency  range  of  200  to  300  Hz  and  a  frequency  step  size  of  0.5  Hz. 


9  Conclusions 


The  Frequency  Window  Reduction  methodology  has  been  implemented  into  a  finite  ele¬ 
ment  environment.  The  prototype  algorithm,  developed  in  F90,  has  been  tested  using  simple 
FEM  models  generated  by  the  ABAQUS™  FEM  code.  The  validity  and  accuracy  of  the  method 
has  been  demonstrated,  with  excellent  results  as  compare  to  full  DOF  model  solutions  even  in 
the  presence  of  numerous  eigenvalues.  A  number  of  the  analytic  features  of  the  FWR  methodol¬ 
ogy  have  been  demonstrated  and  indicate  the  utility  of  these  features  as  analysis  tools  to  extend 
and  enhance  frequency  analysis  as  applied  to  structural  acoustics  problems.  Aspects  of  the  FWR 
methodology  not  only  have  the  potential  to  improve  frequency  analyses  but  also  to  greatly  re¬ 
duce  the  costs  associated  with  these  types  of  analyses. 

The  FWR  codes  having  been  developed  in  standard  F90  can  easily  be  ported  to  other 
computational  platforms  and  may  also  be  extended  into  a  parallel-computing  envirorunent.  The 
most  general  formulation  of  the  FWR  method  includes  complex  frequencies,  complex  frequency 
interpolation  windows  tmd  damping.  The  matrices  (mass  and  stiffness)  which  may  be  obtained 
from  third  party  FEM  codes  such  as  ABAQUS  allow  existing  FEM  codes  to  generate  the  models 
directly  for  the  FWR  method. 

Re-ordering  according  to  a  user-defined  input  “master”  node  list  easily  allows  variable 
degrees  of  reduction  of  models.  This  offers  a  means  by  which  the  optimum  reduction  for  a  given 
model  can  be  determined  and  preliminary  analyses  to  be  performed  in  order  to  characterize  the 
problem.  The  frequency  interpolation  window  location  and  size  can  be  selected  as  well  as  the 
variable  placement  of  internal  interpolation  points  within  the  selected  window  in  order  to  inves¬ 
tigate  and  then  narrow  regions  of  particular  interest  for  the  resultant  system  response.  These 
along  with  full  frequency  window  sweeps  across  the  interpolation  window  or  partial  frequency 
window  sweeps  within  interpolation  window  can  further  refine  the  analysis.  The  choice  of  inter¬ 
polation  methods  can  be  used  to  examine  the  system  response  at  appropriate  levels  of  accuracy. 
The  Lagrangian  interpolation  method  is  generally  less  computational  work  with  a  slightly  re¬ 
duced  accuracy.  While  the  Hermitian  interpolation  method  is  more  computational  work  it  shows 
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highly  accurate  results.  The  eigenvalues  of  the  full  DOF  system  can  also  be  approximated  from 
the  reduced  DOF  system  response,  thus  allowing  the  analyses  to  be  done  without  necessarily  re¬ 
quiring  an  eigenvalue  analysis  of  the  full  DOF  model.  The  combination  of  multiple  interpolation 
windows  into  a  larger  analysis  window  shows  that  the  system  response  of  individual  smaller 
windows  can  be  used  to  build  a  larger  frequency  range  response  tailored  to  the  appropriate  levels 
of  accuracy  from  smaller  windows. 

The  re-use  of  calculated  matrices  is  an  important  aspect  of  the  FWR  method,  which  al¬ 
lows  multiple  analyses  to  be  performed  without  having  to  recalculate  the  necessary  matrices  each 
time.  For  example,  given  un-ordered  mass,  and  stiffness  (and  damping)  matrices,  e.g.  from  a 
third  party  FEM  code,  it  is  possible  to  examine  different  levels  of  DOF  reduction  by  changing 
the  input  “master”  node  list.  Given  this  re-ordered  set  of  matrices,  a  selected  interpolation  win¬ 
dow  is  used  to  calculate  the  constructor  matrices,  which  are  then  valid  for  any  frequency  range 
within  the  defined  interpolation  window.  That  is  the  constructor  matrices  do  not  need  to  be  recal¬ 
culated  as  long  as  the  frequency  sweep  is  within  the  window.  Multiple  frequency  sweep  analyses 
at  different  resolutions  can  then  easily  be  performed.  Further  for  a  set  of  constructor  matrices  de¬ 
fined  for  an  interpolation  window,  the  reduced  DOF  system  matrices  may  be  used  to  perform 
analyses  for  multiple  loading  and  boundary  conditions  without  having  to  recalculate  the  con¬ 
structor  matrices.  This  would  allow  for  a  given  model  variable  load  analyses  to  be  done  quickly 
without  involving  a  major  effort  in  recalculation. 

The  varied  aspects  of  the  FWR  mefiiod  clearly  have  the  potential  to  significantly  impact 
the  way  frequency  analyses  are  performed  on  large  complex  models  as  applied  to  structural 
acoustics.  The  development  and  implementation  of  this  method  as  computational  tools  now  offer 
a  more  effective  and  efficient  computational  environment  in  which  frequency  analyses  can  be 
performed.  The  application  of  these  tools  should  produce  a  major  reduction  of  computational 
costs  while  offering  enhanced  computational  frequency  analyses  of  large  complex  DOF  models. 
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10  Future  Work 


The  Frequency  Window  Reduction,  FWR,  prototype  algorithm  will  continue  to  be  devel¬ 
oped  and  refined  as  a  computational  tool.  The  testing  of  larger,  more  complex  2D  and  3D  struc¬ 
tural  acoustic  models  will  be  undertaken.  This  will  require  the  FWR  codes  to  be  ported  to  a 
UNIX  environment.  With  the  implementation  of  the  FWR  codes  into  a  UNIX  environment,  it  is 
hoped  that  other  research  and/or  industrial  organizations  will  assist  in  the  further  testing  and  de¬ 
velopment  of  this  methodology. 

The  further  investigation  of  the  interpolation  functional  behavior  is  needed.  The  under¬ 
standing  of  the  interdependencies  of  the  degree  of  model  reduction,  the  interpolation  window  lo¬ 
cation  and  size,  the  frequency  sweep  step  size  and  the  choice  of  the  “master”  node  points  is  nec¬ 
essary  in  order  to  refine  the  FWR  method  and  its  overall  efficiency  and  accuracy.  The  usage  and 
utility  of  the  “slave-master”  interpolation  matrix,  a,  as  well  as  more  arbitrary  loading  and  con¬ 
straint  conditions  will  be  investigated  as  applied  to  larger  FEM  models.  The  matrix  properties  of 
the  constructor  and  P  matrices  need  to  be  characterized  in  order  to  develop  and  implement  more 

efficient  and  optimize  storage  and  computational  schemes.  Techniques  to  obtain  the  p  matrix  di¬ 
rectly  and  implement  the  interpolation  of  this  matrix  by  similar  methods  as  for  G”‘will  be  ex¬ 
plored.  This  would  allow  more  arbitrary  loading  and  constraint  conditions  on  the  “slave”  DOFs 
and  also  lead  to  reconstruction  of  the  full  L"‘  system  operator  from  the  reduced  system  operator 


The  utility  and  implementation  of  this  methodology  will  be  investigated  for  other  fre¬ 
quency  dependent  applications  such  as  electromagnetic  FEM  modeling  and  inverse  scattering 
problems.  Finally,  the  direct  incorporation  of  the  resonance  or  modal  response  into  the  compu¬ 
tational  methodology  will  be  undertaken  and  tested  on  FEM  models. 
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12  Appendices 

Appendix  A  -  Example  of  Master  and  Slave  Degrees  of  Freedom 

As  described  above,  the  “master”  DOFs  are  chosen  as  DOFs  which  have  applied  loads  or 
fixed  boundary  conditions.  In  addition  to  loaded  and  constrained  DOFs,  the  master  DOFs  may  be 
chosen  as  important  place-holders  represented  in  the  models;  the  interfaces  or  boundaries  be¬ 
tween  element  groups  of  differing  material  properties  or  element  grouping  midpoints.  The 
“slave”  degree  of  freedom  refers  to  the  dependent  DOFs  or  those  DOFs,  which  are  to  be  con¬ 
densed  or  reduced  out  of  the  full  FEM  model.  These  are  DOFs  which  typically  do  not  have  ap¬ 
plied  loads  or  constraints  and  do  not  contribute  to  the  system  response  or  are  not  of  interest. 

A  simple  DOF  example  base  on  an  FEM  model  is  illustrated  in  the  figure  below 


304S  Steel  Poly: 
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•  Master  Node  Points 

In  this  FEM  model  it  is  assumed  that  of  the  46  node  points  there  are  14  “master”  node 
points  and  so  there  would  be  32  “slave”  node  points.  Here  the  “master”  node  points  are  selected 
as:  node  points  which  have  the  applied  load,  the  fixed  boundary  condition  ,  the  boundaries  be¬ 
tween  element  groups  of  differing  material  properties  and  the  element  grouping  midpoints.  The 
original  FEM  model  nodal  numbering  can  then  be  rewritten  via  a  reordering  look-up  table  as 
listed  below.  This  table  now  allows  conversion  between  the  original  node  numbering  as  is  typical 
in  FEM  models  and  the  “master”  node  numbering  used  in  the  CMR  method. 
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Node  Numbering  Look-Up  Table 


I 

n 

HI 

Node 

“Master”  Node  Numbering 

Original  Node  Numbering 

Number 

In  Original  Index 

In  “Master”  Index 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

1 

2 

15 

16 

17 

18 

19 

20 

5 

6 

21 

22 

23 

24 

15 

3 

25 

16 

4 

26 

17 

5 

7 

18 

6 

8 

19 

7 

27 

20 

8 

28 

21 

11 

29 

22 

12 

30 

23 

13 

9 

24 

14 

10 

25 

15 

31 

26 

16 

32 

27 

19 

33 

28 

20 

34 

29 

21 

11 

30 

22 

12 

31 

25 

35 

32 

26 

36 

33 

27 

37 

34 

28 

38 

35 

31 

39 

36 

32 

40 

37 

33 

13 

38 

34 

14 

39 

35 

41 

40 

36 

42 

41 

39 

43 

42 

40 

44 

43 

41 

45 

44 

42 

46 

45 

43 

3 

46 

44 

4 
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For  the  above  example  the  “master”  node  numbering  is  input  as  a  list  of  14  node  numbers 
from  the  FEM  model  which  represent  the  selected  nodes  which  are  to  be  the  “master”  nodes.  ( 
Shaded  area  in  column  II  of  the  table  )  The  remainder  of  the  list  (Column  II,  #s  15  to  46)  which 
is  automatically  generated  represents  the  “slave”  nodes,  which  will  be  condensed  out  of  the 
model.  The  “master”  or  reordered  node  numbering  correspondence  to  original  or  FEM  model 
node  numbering  which  is  also  automatically  generated  is  listed  in  column  III.  The  interpretation 
of  this  column  stated  in  words  is  the  following...  for  example  node  number  3  in  the  original 
FEM  model  (column  I)  is  found  via  column  III  to  be  positioned  as  node  number  15  in  the  reor¬ 
dered  numbering  (column  H).  This  look-up  table  can  now  be  used  to  either  directly  assembles  the 
global  matrices  in  the  reordered  numbering  scheme  or  to  reorder  existing  global  matrices. 

Appendix  B  -  Frequency  Coordinate  Transformation 

2D  Case 

If  the  frequency,  co  ,  is  assumed  to  be  a  complex  variable,  m  =  O),  +  itOj  where  the 
subscript  “1”  refers  to  the  “real”  part  of  co  and  “2”  refers  to  the  “imaginary”  part  of  co,  then 
q  is  defined  as  a  complex  variable,  q  =  q^  +  iq2  ■  It  is  then  possible  to  define  the  variable 
transformation  q  =  io)  and  such  that  q  =  —0)2  +  ico^  ■  Hence,  q^  =—0)2,  i.e.  q^  is  the 
“real”  part  of  q,  and  ^2  =  ^2  is  the  “imaginary”  part  of  ^  . 

A  rectangular  window  in  ©-space  (frequency  window)  can  be  defined  by  the  “real”  and 

“imaginary”  origin  O)]  ,  i(02  and  the  window  width,  as  defined  along  the  “real”  axis  and 
height,  as  defined  along  the  “imaginary”  axis.  Hence,  the  four  vertices  of  the  rectangular 
window,  cOf ,  are 

dJj  =  Q)^  +  ia>2 

©2  =  (®i+W„)  + 

©3  =  (®, +W<^)+  i(©2+/i„) 

©4  =  ©1  -t-  1(0)2 +h^) 

Note:  ©,,  ©2  and  the  window  height  h^,  and  width  w^,  are  all  real  values  in  ©-space . 
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Via  the  coordinate  transformation  q  =  ifi) ,  a  frequency  window  may  now  be  referred  to 
in  q-space,  with  an  origin  q^,  q^  and  a  window  height,  h^,  and  width,  w^.  The  window  coordi¬ 
nate  origin  is  then  q  =  i(n),  -I-  i(Oj)  =  -(0^+  ,  and  so  the  window  vertices  are 

q^  =  =-6)2+ 

^2  =  *6)4  ^-((02+K)  +  itOi 

^3  =  1^3  =-(t02+^<a)  + 

^4=  iWj  =-t»2  + 

Note:  ^2  the  window  height,  h^,  and  width,  ,  are  all  real  values  in  q-space. 

Given  that  =  -m2and^2  =  ft),  and  also  =  h^mdh^  =  the  window  vertices  are 


^1=^1+  i  (I2 

h  =  +  *  ^2 

q,  =  (^i+wj  +  i  (q2+\) 
q,=  qi  +  i  (^2+\) 


The  coordinate  transformation  q  =  ico  ,  may  be  simply  thought  of  as  a  rotation  of  the  fre¬ 
quency  axis  by  90°. 


Imaginary  O) 


Imaginary  q 
92^  ®i 


q,  =  -  (O2 

◄ - 

Realq 


ID  Case 

For  the  ID  case,  the  interpolation  points  of  a  frequency  range  or  line, 
Wj (ft), ,(02’^e>’K)  ~  defined  as  the  line  end  points 
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(Do  =  6), 

+  i<D2 

either 

m,  =  (<0, 

+  w^)+  i(02 

if  K  =  o 

or 

m,  =  m, 

+  i{(02  +  K) 

if  w<j,=  0 

Imaginary  co 

I  T' 

• - •  K=0 

a 

•  % 

- ►  Real  CO 

As  above  for  the  2D  case  via  the  coordinate  transformation  q  =  ico  ,  a  frequency  line 
may  now  be  referred  to  in  q-space,  with  an  origin  q^,  q^  and  a  window  height,  h^,  and  width, 
w^.  The  line  coordinate  end  points  in  q-space  are  then  found  as 

q^  =  iuJo  =  -CO2 

q^  =  im^  =-{(02 +K)  +  *  if>»^co=0 

q,=  =-©2+  i(cui+w^)  ifK  =  0 

Appendix  C  -  Reduced  Mass,  Stiffness  and  Damping  Matrices 

A  frequency  m  =  m,  +  iOj  is  transformed  to  the  generalize  variable  q  =  q^  +  iq2 
hy  q  =  i<5) .  Then  for  a  given  value  of  ^  ^2 ) '''ithin  the  interpolation  window,  the 

matrix  may  be  founds  as 

IS(q)  =  G-'(«)  H{q) 

Where  Hiq)  =  a{q)  (^^M„  +  ^D„  +  S„)  -  (^^M2i  +  ^D2,  +  S21). 
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The  G“'(^)  matrix  for  Hermitian  Interpolation  for  a  given  value  of  g  within 

the  window  having  previously  obtained  the  for  the  rectangular  window  at  the  points 

is 

G-'(q)  =  C„(«„4)  W(q,) 

Alternatively  the  value  of  G~^{q)  can  be  interpolated  by  the  Shape  Functions  evaluated  at 
q  =  /(4,4)  from  the  isoparametric  coordinates  {r,s)  and  the  matrices  previously 

calculated  at  each  of  the  “element”  or  window  interpolation  points  by 

=  'Zhi{r,s)  G;\q) 

Finally  reduced  system  operator,  Lo(^)  is 

L,{q)  =  p{q) 

=  +  ^1^11  +  ^u)  +  Sjj)  P{q) 

From  this  equation,  the  reduced  mass,  stiffness  and  damping  matrices  may  be  obtained  as 
Mo(^)  =  M„  +  M,,  P{q) 

So(^)  =  s„  +  s„  m 

Do(^)  =  D„  +  D„  p{q) 

Hence  in  order  to  calculate  the  reduced  mass,  stiffness  and  damping  matrices  the  matrix 
P{q)rmst  be  obtained  at  each  specific  value  of  the  frequency  (or  q). 
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